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Appendix. Mitrovica et al. [1994a] and [1994b] 



I. Introduction 


This project involves obtaining GPS measurements in Scandinavia, and using the 
measurements to estimate the viscosity profile of the Earth’s mantle and to correct 
tide-gauge measurements for the rebound effect. Below, we report on several aspects 
of this project . 

II. GPS Measurements 

The DSGS was not scheduled to be reoccupied with DOSE receivers during the 
report period. The permanent network set up by Onsala Space Observatory continues 
to operate, and the data are being evaluated. 

III. Theoretical Advances 


An important technical advance we intend for this project is to use the full three- 
dimensional site velocity information for inferring geophysical parameters. During 
the report period, two papers [Mitrovica et al. 1994a,b] have been been accepted for 
publication in the Journal of Geophysical Research and will be published in April. 


Reprints of these papers are contained in the Appendix. 
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A spectral formalism for computing three-dimensional 
deformations due to surface loads 

I. Theory 

J. X. Mitrovica 

Department of Physics, University of Toronto, Toronto, Ontario, Canada 

J. L. Davis and I. I. Shapiro 

Harvard-Smithsonian Center for Astrophysics, Cambridge, Massachusetts 


Abstract. We outline a complete spectral formalism for computing high spatial 
resolution three-dimensional deformations arising from the surface mass loading 
of a spherically symmetric planet. The main advantages of the formalism are 
that all surface mass loads are always described using a consistent mathematical 
representation (expansions using spherical harmonic basis functions) and that 
calculations of deformation fields for various spatial resolutions can be performed by 
simply altering the spherical harmonic degree truncation level of the procedure. The 
latter may be important when incorporating improved observational constraints on a 
particular surface mass load, when considering potential errors in the computed field 
associated with mass loading having a spatial scale unresolved by the observational 
constraints, or when treating a number of global surface mass loads constrained 
with different spatial resolutions. These advantages do not extend to traditional 
“Green’s function” approaches which involve surface element discretizations of the 
global mass loads. In the practical application of the Green’s function approach the 
discretization is subjective and dependent on the particular surface mass load being 
considered. Furthermore, treatment of mass loads with higher spatial resolutions 
can require tedious rediscretization of the surface elements. Another advantage 
of the spectra] formalism, over the Green’s function approach, is that a posteriori 
analyses of the computed deformation fields, such as degree correlations, power 
spectra, and filtering analyses, are easily performed. In developing the spectral 
formalism, we consider specific cases where the Earth’s mantle is assumed to 
respond as an elastic, slightly anelastic, or linear viscoelastic medium. In the case 
of an elastic or slightly anelastic mantle rheology the spectral response equations 
incorporate frequency dependent Love numbers. The formalism can therefore be 
used, for example, to compute the potentially resonant deformational response 
associated with the free core nutation and Chandler wobble eigenfunctions. For 
completeness, the spectral response equations include both body forces, as arise 
from the gravitational attraction of the Sun and the Moon, and surface mass 
loads. In either case, and for both elastic and anelastic mantle rheologies, we 
outline a pseudo-spectral technique for computing the ocean adjustment associated 
with the total gravitational perturbation induced by the external forcing. Three- 
dimensional deformations computed using the usual Love number approach are 
generally referenced to an origin located at the center of mass of the undeformed 
planet. We derive a spectral technique for transforming the results to an origin 
located at the center of mass of the deformed planet. 


Copyright 1994 by the American Geophysical Unions 
Paper number 93JB03128. 
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Introduction 


The Earth is subject to a number of external forcings 
having characteristic time scales which vary over many 
orders of magnitude. The gravitational potential due to 
the Sun and the Moon gives rise to both body tides and 
surface mass load variations associated with ocean tides 
which have significant power at periods ranging from less 
than a day to several decades [e.g., Cartwright and Tayler, 
1971], Time series of surface pressure measurements in- 
dicate comparably large ranges in atmospheric loading 
time scales [e.g., Trenberth, 1981; Mtiller and Ziim, 1983; 
Rabbel and Zschau, 1985]. Furthermore, the mass bal- 
ance of continent-based ice sheets fluctuates with time 
scales ranging up to the period of approximately 10 s years 
characteristic of the late Pleistocene glacial cycles [e.g., 
Shackle ton, 1967; Broecker and Van Donk, 1970], though 
small land-based ice sheets and glaciers are known to 
be presently melting over significantly shorter timescales 
[ Meier, 1984]. 

The prediction of solid Earth deformations associated 
with these surface mass loads is a classic problem in geo- 
physics. (See Farrell [1972] for a detailed discussion.) The 
advent and dramatic improvements in the accuracy of 
space geodetic measurement techniques, including very 
long baseline interferometry (VLBI), surveying using the 
Global Positioning System (GPS) and satellite and lunar 
laser ranging (SLR and LLR), have revived considerable 
interest in such predictions. Indeed, analyses of the ef- 
fect of atmospheric loading [e.g., van Dam and Wahr, 1987; 
van Dam, 1991], ocean loading [e.g., Schemeck, 1991], and 
late Pleistocene glacial loading [e.g., Tushingham, 1991] on 
baselines determined using space geodetic measurement 
techniques have all been recently considered. 

Solid Earth deformations associated with surface mass 
loads have been computed using an approach developed 
from the work of Longman [1962, 1963] and Farrell [1972]. 
Farrell [1972], following Longman [1962, 1963], computed 
Green’s functions for the radial and tangential displace- 
ments associated with the impulse loading of a spherically 
symmetric planet having an elastic mantle rheology. The 
response of the Earth model to an arbitrary surface mass 
load can then be computed via a surface convolution of 
the mass load with the appropriate Green’s function. In 
the practical application of Farrell’s [1972] Green’s func- 
tion approach the surface mass load is commonly dis- 
cretized using a number of surface elements [Jentzcsh, 
1985]. The response of the Earth model as a function 
of position relative to a set of such surface elements with 
varying geometries is precomputed. The two-dimensional 
convolution required to obtain the response to Em arbi- 
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tr&ry surface mass toad then reduces to a weighted sum of 
a relatively small number of “surface element” responses. 
The discretization procedure clearly reduces the effort in- 
volved in generating a number of different deformation 
fields. Nevertheless, it is also motivated by the finite 
spatial resolution with which the surface mass loads are 
known. As an example, van Dam and Wahr [1987] dis- 
cretized their atmospheric load into a 2.5° x 2.5° grid 
system to be consistent with the resolution of the meteo- 
rological data sets available to them. 

There have been recent attempts to compute Green’s 
functions for Earth models of increasing complexity. For 
example, a number of analyses have incorporated the ef- 
fect of mantle anelasticity on the impulse response [Wang, 
1986; Schemeck, 1991; Pagiatakis, 1990]. Furthermore, Pa- 
giatakis [1990] has also considered the influence of man- 
tle anisotropy. In a summary of these results, Schemeck 
[1991] has shown that discrepancies between the Green’s 
functions derived by Schemeck [1991] and Pagiatakis [1990] 
and those of Farrell [1972] are generally less than 2%. 

The Green’s function approach has also been extended 
to consider the (deformational and sea level) response of 
an Earth model with a linear viscoelastic mantle rheology 
to an applied surface mass load [Peltier, 1974]. The most 
common application is the analysis of the response of the 
Earth to the last major deglaciation event of the current 
ice age [e.g., Peltier and Andrews, 1976; Clark et aL, 1978; 
Wu and Peltier, 1983; Peltier, 1985; Tushingham, 1991; James 
and Lambert, 1993]. In the case of a viscoelastic mantle 
rheology, the impulse response Green’s functions are de- 
pendent on time, as well as on the angular distance from 
the point load. As a consequence, a calculation of the 
response of an Earth model to various surface load ele- 
ments becomes substantially more time consuming. In 
this regard, the recent analysis of Tushingham and Peltier 
[1991] has adopted a circular disk load discretization of 
the late Pleistocene ice sheets in which the smallest disk 
diameter is 1° . 

In its practical application the main disadvantages of 
the Green’s function approach are that the surface ele- 
ments are arbitrarily chosen (and are often dependent, in 
the literature, on the particular global mass load being 
considered) and that changes in the spatial resolution of 
the surface mass load often require a tedious rediscretiza- 
tion of the load geometry. (Compare, for example, the 
discretizations of Wu and Peltier [1983] and Tushingham 
and Peltier [1991].) It is, as a consequence, difficult to 
consider the effect on the predicted deformation field of 
discretization errors associated with the representation of 
the surface mass load. Furthermore, the a priori construc- 
tion of the Green’s function generally involves a sum to 
very high harmonic degree ( Farrell [1972] truncated the 
sum at degree 10,000) which is not necessarily consistent 
with the observationally constrained resolution available 
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regarding the surface mass load. In the case of an elas- 
tic mantle rheology this fact is of little consequence since 
the summation is performed only once. However, in the 
viscoelastic case, where the radial viscosity profile is gen- 
erally a free parameter, Green’s functions (which depend 
on this profile) may be computed many times. Hence a 
method which incorporates the impulse response of the 
Earth model to a spherical harmonic degree consistent 
with the resolution of the surface load is advantageous. 

In this paper we outline a spectral, i.e., spherical har- 
monic, formalism for computing three-dimensional defor- 
mations arising from surface mass loads. Within the for- 
malism, all surface mass loads are described using a con- 
sistent mathematical representation based upon spherical 
harmonic decompositions and the surface convolutions re- 
quired to compute the radial and tangential surface dis- 
placements associated with the load can be performed an- 
alytically. In this case, the resulting analytic expressions 
truncate at the spherical harmonic degree level chosen for 
the surface mass load (which, in turn, will be governed 
by the observational constraints on that mass load). Con- 
sideration of different spatial resolutions will require no 
more preliminary effort than a change in this truncation 
level. Another obvious advantage of a spectral formalism 
is that a posteriori analyses of the computed deformar 
tion fields, such as the determination of power spectra 
and degree correlations or the filtering of specific spatial 
wavelengths, are trivially performed. 

As Farrell (1972] has pointed out in regard to the 
surface-mass-load problem, the primary disadvantage of 
spectral techniques is that their spatial resolution is in- 
dependent of position. As a consequence, incorporating 
highly resolved observational constraints at one location 
on the Earth’s surface forces better resolution than may 
be necessary from observational constraints from other 
regions. The computational power and availability of 
modem work stations, however, make this argument less 
of a practical limitation than it was 20 years ago. For 
example, we have implemented the algorithms presented 
below on a Hewlett-Packard model 730 work station (em- 
ploying RISC architecture) with 32 Mbytes of RAM and 
found that with no special techniques we are able rapidly 
to obtain spherical harmonic truncation levels of degree 
2048 when considering the deformation field induced by 
a specific surface mass load. This resolution corresponds 
to a spatial resolution of 0.1° or 10 km, superior to the 
resolution adopted in analyses of the deformation of the 
Earth based on the surface element Green’s function ap- 
proach. (The spatial resolutions described above in the 
context of recent analyses of this type are an indication 
of this disparity.) In this regard, the spatial resolutions 
we have obtained are certainly better than those which 
characterize current observational constraints on the var- 
ious global surface mass loads. These results suggest that 
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it is time to reconsider the utility of spectral techniques 
in analyses of Earth deformations by surface loads; it is 
intended that this study will provide the algorithms nec- 
essary to do so. 

To our knowledge, the analysis of Pertsev [1966] was 
the first to suggest using a spectral approach to com- 
pute the response of an Earth model whose Love numbers 
have been calculated. Pertsev [1966] was specifically con- 
cerned with variations in gravity associated with semidi- 
urnal ocean tides on an elastic Earth. In addition, Peltier 
[1976] derived, but did not apply, a spectral formalism 
for computing radial deformations on a linear viscoelas- 
tic Earth. Spectral formalisms have also been derived 
to compute the redistribution of ocean mass associated 
with the melting of late Pleistocene ice sheets on a vis- 
coelastic Earth [Nakada and Lambeck, 1987; Mitrovica and 
Peltier, 1991a], The Nakada and Lambeck [1987] approach 
represents an approximate solution to the governing sea 
level equation, while the spectral and pseudo-spectral al- 
gorithms described by Mitrovica and Peltier [1991a] yield a 
gravitationally self-consistent solution of the equilibrium 
ocean adjustment. Finally, Dickman [1989] has outlined 
a spherical harmonic approach to predicting long-period 
ocean tides which incorporates a number of important 
physical effects (ocean bottom friction, etc.). Regardless 
of the assumed mantle rheology, horizontal deformations 
induced by surface mass loads have not, to date, been 
treated using a spectral approach. 

In deriving a spectral formalism for computing three- 
dimensional deformations, we will consider Earth models 
characterized by elastic, slightly anelastic, and viscoelas- 
tic mantle responses. In the case of elastic and anelastic 
mantle rheologies, we incorporate frequency dependent 
Love numbers in the derivation of the spectral response 
equations. As a consequence, the formalism can be used, 
for example, to predict the resonant deformational re- 
sponse associated with mass load variations having fre- 
quencies near the eigenfrequencies of the free core nutar 
tion and Chandler wobble normal modes of the planet 
[Wahr, 1981a, b]. The spectral response equations require, 
as input, the spherical harmonic coefficients of the surface 
mass load and the Love numbers appropriate to the rhe- 
ological model being considered. The reader is referred 
to the discussions of Longman [1962, 1963] and Farrell 
[1972], Wahr and Bergen [1986], and Peltier [1974, 1985], 
for the treatment of Earth models with, respectively, elas- 
tic, anelastic, and linear viscoelastic mantle rheologies. 

The Love numbers are generally defined such that dis- 
placements computed using them are referenced to an 
origin located at the center of mass of the undeformed 
planet. (See Farrell [1972] for a detailed discussion of 
this point.) As a consequence, the computed positions 
will include a component associated with a pure translar 
tion of the center of mass of the deformed planet (which 
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does not include the surface mass load) relative to its 
original position. To complete the present study, we pro- 
vide a straightforward spectral approach for transforming 
the computed positions to a reference frame whose origin 
is located at the center of mass of the deformed planet. 
The latter positions are more accurately interpreted as 
“deformations.” 

For the case of elastic and anelastic mantle rheologies 
we provide a pseudo-spectral technique for computing the 
ocean adjustment associated with the gravitational per- 
turbation induced by the surface mass loading. The tech- 
nique outputs the spherical harmonic coefficients of the 
ocean adjustment, as do the algorithms of Nakada and 
Lambeck [1987] and Mitrovica and Peltier [1991a] for the 
viscoelastic case, as well as that of Die km an [1989] for 
his ocean tide theory. As a consequence, the spectral re- 
sponse equations derived herein are well suited for com- 
puting the three-dimensional crustal deformations asso- 
ciated with these mass redistributions. 

In a companion article [Mitrovica et al., this issue] we 
considered predictions of present-day three-dimension- 
al crustal deformation rates associated with the melt- 
ing of the late Pleistocene ice sheets. We will use the 
spectral response equations derived herein to predict the 
three-dimensional deformations which arise from a real- 
istic space-time deglaciation chronology and gravitation- 
ally self-consistent ocean mass adjustments. 

A Love Number Formulation of the Earth’s 
Response to External Loads' 

In this and following sections we develop a complete 
spectral formalism for computing the response of a spher- 
ically symmetric, self-gravitating planet to an arbitrary 
external load. We focus mainly on radial and tangen- 
tial surface deformations. Nevertheless, the theory can 
also be used with the appropriate Green’s function (see 
below) to compute anomalies in the gravitational field 
of the planet. As discussed in the introduction, we con- 
sider elastic, anelastic, and linear viscoelastic mantle rhe- 
ologies and in each case incorporate gravitationally self- 
consistent ocean loading components in the formulation. 
The spectral formalism is based on the usual Love num- 
ber formulation of the external loading problem which we 
review, briefly, below. 

Farrell [1972], following Longman [1962, 1963], outlined 
a numerical procedure for determining the tidal and sur- 
face load Love numbers of a spherically symmetric, self- 
gravitating Earth having an elastic mantle and a fluid 
core. These numbers, originally defined by Love [1909] 
(and also by Shida and Matsuyama [1912]), govern the re- 
sponse of the planetary model to a gravitational potential 
and surface mass load forcing, respectively. In consider- 
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ing the Earth’s tidal response, let us begin by assuming an 
Earth model which experiences no frictional dissipation 
(i.e., we assume no anelastic component of the response). 
In this case the tidal Love numbers may be written in the 
form 

hj = hJ' E (1) 

If = lJ’ E ( 2 ) 

kj = kJ' E (3) 

where the superscript E denotes a coefficient determined 
from a model with a purely elastic mantle rheology and l 
is the spherical harmonic degree. The Love numbers hj, 
ij, and kj represent coefficients in harmonic expansions 
of the radial and tangential displacement and a poten- 
tial perturbation due to the planetary mass redistribu- 
tion (i.e., body tide), respectively, induced by the tidal 
potential. 

Let us denote the tidal potential on the Earth’s surface 
at a geographic location of colatitude 6, east longitude 
rp, and time t, as W (6, ip, t). Practical representations 
of W(0,ip y t) have traditionally been based on the de- 
termination of time dependent coefficients in a spherical 
harmonic expansion of the tidal potential [e.g., Cartwright 
and Tayler, 1971; Cartwright and Edden, 1973]. We can thus 
write 


OO / 

W(6,rP,t) = X! E "WO • ( 4 ) 

/=2 m——t 

The Y( m (0, ip) are surface spherical harmonics normalized 
(in our formalism) such that 

JJ Y S>T n '( 6 >'P) Y ern(0,ip) sine dOdrp = 
n 

4n6 u >6 mm i , (5) 

where JJ n is used to indicate integration over the entire 
solid angle and % is the Kronecker delta. 

Using (l)-(3) and (4), the response of the Earth to the 
tidal potential may be written as [Melchior, 1983] 

1 00 £ 

U T (0, xp t t) = - £ 2 hJ’ E WWO Y tm (0, xp) (6) 

9 t=2 m=-l 

? T (0,tM) = ;E E l t’ B V Y im(0,tp) (7) 

9 1=2 m=-l 

and 


4> T {0,xp,t) 


OO l 


E E k ? E 


Y tm (0,rP) , 


( 8 ) 
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where U 7 and V T are the radial and tangential body 
tide displacements, <p 7 is the potential perturbation due 
to the planetary mass redistribution, and g is the surface 
gravitational acceleration. The symbol V denotes the 
two-dimensional gradient operator: 


V = 


- d 

e — i-t b 
do w 


_i d_ 

sin 0 dtp ' 


(9) 


The terms VYi m (0, ip) will play an important role in this 
study, and in the appendix we provide expressions for 
them in terms of a set of Ye m (0,ip). 

Equations (6)— (8) may be used to compute the bar 
sic deformational and gravitational response of the solid 
Earth to the tidal potential forcing. The equations can 
be combined to generate expressions for other responses 
of geophysical importance. For example, in considering 
the gravitationally self-consistent redistribution of ocean 
mass induced by the tidal potential (as we do below), 
it is useful to have an expression for the gravitational 
potential perturbation at the deformed solid surface in- 
duced by both the tidal potential and the resulting body 
tides. The perturbation, which we denote by $ 7 (0,ip,t), 
is [Longman, 1963] 


* T (0,ip,t) = 


OO £ r 

££ l + 

*= 2 m=— / L J 


W tm (t) Y tm ($,4>) , (10) 


where the term equal to unity reflects the direct effect of 
the tidal potential and the term h 7 ' E arises from trans- 
forming from the undeformed to the deformed solid sur- 
face. 

Let us turn now to the surface load Love numbers. In 
computing the response of the Earth to variations in sur- 
face mass loads we will consider examples where both 
elastic and linear viscoelastic mantle rheologies are rele- 
vant. As a consequence, we require an extension of the 
elastic Love number forms (1)— (3) to the general linear 
viscoelastic case. Such an extension has been provided 
by Peltier [1974, 1976], who invoked the correspondence 
principle [Biot, 1954] to derive so-called viscoelastic sur- 
face load Love numbers. These dimensionless numbers 
define the response of the Earth model to an applied point 
impulse load. In the time domain these Love numbers, 
which we denote as hf, l£, and k E , have the following 
normal mode form [Peltier, 1976] 

K(t) 

ti(*) = hf^t) + £) rlj^expf-^*) (11) 

*=i 

K(t) 

#(«) = t' E m + £ r2^exp (sit) (12) 

*=i 
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and 


K(t) 

k{(t) = k E ’ E 6(t) + Y, r3 k‘ exp(-sit) . (13) 

fc=i 

The first term on the right-hand sides of each of (11)- 
(13) represents the instantaneous elastic response to an 
applied point mass impulse load. The second term is 
the nonelastic response which is represented by a super- 
position of K(£) modes of pure exponential decay. The 
formalism outlined by Peltier [1985] can be used to com- 
pute the normal mode amplitudes rl£’*, r2 E,t , and r 3%’* 
and the inverse decay times (eigenfrequencies) s%. 

The hf. If, and kf surface load Love numbers repre- 
sent the coefficients of degree l in Legendre polynomial 
expansions of Green’s functions for the radial and tan- 
gential displacement and the potential perturbation due 
to mass redistribution of the planet, respectively. The 
Green’s functions can be written as 

OO 

g u ( 7 ’«) = t t ) p *( cos 7) ( 14 ) 

e 1=0 
OO 

Gy(7,i) = Jj- 'Y'ttW gZ p t( cos 7) 7 ( 15 ) 

e /=0 7 

and 

OO 

G$(%0 = W e J^ 0 ^ {t) p ‘( cos ^' ( 16 > 

respectively, where o is the mean radius of the Earth and 
M e its mass, 7 is a unit vector tangent to the surface (i.e., 
horizontal) in the direction of the great circle extending 
from the load point (O', rp') to the observation point ( 0 , rp), 
and 7 is the angular distance between these two points, 
given by 

cos 7 = cos 0 cos O' -f sin 0 sin O' cos (rp — rp') . (17) 

The Green’s functions (14)— (16) can be combined to 
generate Green’s functions for a number of geophysical 
observables. Mitrovica and Peltier [1989] have derived 
Green’s functions for both the free air gravity anomaly 
and the geoid anomaly. (See also Mitrovica and Peltier 
[1991b].) Longman [1963] determined Green’s functions 
for both the gravity anomaly and the gravitational po- 
tential perturbation at the deformed surface for the case 
of an elastic mantle rheology. In the viscoelastic case the 
gravitational potential perturbation Green’s function 
has the following form [Peltier and Andrews, 1976]: 

g *(7, ‘)-gf:[(i+ k r ,E - h <’ B ) w 

K(t) 

*=1 


( r3 k’ £ ~ rl k’ e ) exp -P/(cos 7 ) (18) 
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where the term equal to unity in the multiplier term of 
the delta function incorporates the direct gravitational 
attraction of the impulse load and the terms in h E,E and 
rl£’* arise from transforming from the undeformed to the 
deformed surface. 

The response of a spherically symmetric, self-gravitat- 
ing, viscoelastic planet to an arbitrary surface mass load 
L(9,ip,t) (with dimensions of mass per unit area) can be 
computed by convolving the load, in the space and time 
dimensions, with the appropriate Green’s function. If we 
denote a general Green’s function as G L (*f,t), then the 
response R is given by 

R{9,tp,t) = 
e 

j jj^w, dt\ (19) 

— oo n 

where 7 is given by (17). The response R(9,tp,t) will 
be a scalar quantity for the case of the Green’s functions 
Gy, G%, and and a vector quantity for the tangential 
displacement Green’s function Gy. 

The Green’s functions (14)— (16) represent extensions to 
the viscoelastic case of the analogous Green’s functions 
described by Longman [1963] and Farrell [1972] for the 
purely elastic case. Indeed, the latter can be derived from 
the former by considering only the elastic component of 
(11)— (13) in the summations of (14)— (16), that is, for the 
special case of a purely elastic response, we can write 

g l v e (i) = ^ E h i’ Ep ‘( c 08 1) ( 20 ) 

* <=o 

G t >E h) = k p’ B r, p ^ (21) 

and 

g * e w = |rfX’ ep * (cos7) - (22) 

In this case the general response (19) would simplify to 

R E (9,rP,t) = JJ o?H9\rP\t)G L ’ E { 7 ) dfl'. (23) 

n 

The now classic method of computing the response 
R E (9,tp,t) discussed by Farrell [1972] consists of a di- 
rect application of (20)-(23). First, the Love numbers 
h E ' E , lj" E , and k E ’ E are computed for a specific “elas- 
tic” Earth model. Second, the summations in (20)-(22) 
are computed to very high spherical harmonic degree by 
making use of asymptotic properties associated with infi- 
nite sums. Finally, the convolution between the Green’s 


10 



function and the specified load L(9,ip,t) is computed in 
the space domain numerically. In this respect, the two- 
dimensional convolution can be discretized by using (20)- 
(22) to compute Green’s functions for finite surface el- 
ements [Jentzsch, 1985]. These elements are commonly 
chosen to be circular disks (in which case “equivalent” 
Green’s functions need to be computed for disks of various 
sizes and at various angular distances from the observa- 
tion point), but Goad [1980], for example, has advocated 
the use of surface ring sectors. Other discretizations are 
also common [e.g., van Dam, 1991]. 

The Farrell [1972] methodology has been adopted by, 
for example, Rabbel and Zschau [1985], van Dam and Wahr 
[1987], and van Dam [1991] in considering deformations 
and gravity changes associated with atmospheric loading. 
The same approach has also been used to consider the 
planetary response to tidal loading [Goad, 1980; Wang, 
1986; Francis and Dehant, 1987; Schemeck, 1991]. 

The majority of numerical predictions of the glacial 
isostatic adjustment of the Earth in consequence of the 
melting of the late Pleistocene ice sheets have proceeded 
via an extension of Farrell's [1972] space-domain convo- 
lution approach to the viscoelastic case, as defined by 
(14)— (19). These studies include the analysis of gravita- 
tionally self-consistent sea level change [Farrell and Clark, 
1976; Peltier et al., 1978; Wu and Peltier, 1983], as well as 
crustal deformations [e.g., Peltier, 1974; James and Mor- 
gan, 1990; Spada et aL, 1992], Indeed, as discussed in 
the introduction, predictions of three-dimensional crustal 
deformations associated with the phenomenon have only 
been computed using such an approach. In the viscoelas- 
tic case, calculation of the point mass Green’s functions 
(or of the surface element Green’s functions) requires sub- 
stantially more effort than in the elastic case, since the 
functions are time dependent. 

In the next section we outline a special approach which 
allows for the analytic determination of the required sur- 
face convolution integrals. In subsequent sections the re- 
sults are applied to consider the response of Earth mod- 
els having elastic, anelastic, and viscoelastic rheologies, 
to arbitrary surface loads. 

Response to an Arbitrary Surface Mass Load 
of a Spherically Symmetric Earth: 

A Spectral Approach 

As in the case of the tidal potential, let us represent 
the surface load L(9, ip, t ) in terms of a spherical harmonic 
expansion of the form 

L(M,t) = X; Limit) Y tm (0,rP) } (24) 

t,m 
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where we have adopted the short-form notation 


OO / 

!>££• 

i,m 1=0 m=—t 

The surface load Green’s functions described in the last 
section are of two forms, which we may denote symboli- 
cally as 

OO 

Gi(7,t)=X>/(*)^( c °s7) (25) 


<3a(7, t) = £ B e (t) £-Pt ( cos 7 ) 7 , (26) 

t=o 

where (14), (16), (18), (20), and (22) are of the form (25), 
whereas the vector Green’s function for the tangential 
displacement (15) and (21) are of the form (26). 

To begin, let us consider the form (26). Using (26) and 
(24) in the response (19) yields 


fl(w)=fy T 

t = 0 


X 



W(t') 


x II Ye ' m ^ d ' , ^ , ^ pe( ' cos ' y ^ 

n 


dt' (27) 


It is straightforward to show that, in general, 

J^/(7)7 = V/(7), (28) 

where V is the two-dimensional gradient operator defined 
in (10), 7 is the arc length defined by (17), 7 is a unit vec- 
tor defined below (16), and /( 7 ) is an arbitrary function. 
Equation (28) indicates that the directional derivative op- 
erator 7 <?/f ?7 acting on /( 7 ) is equivalent to the gradient 
operator V acting on /( 7 ). The gradient operator is, of 
course, independent of the load point coordinate (O', ip'). 
Using (28) in (27) yields 


OO ‘ 

R(0,iP,t) = J^a 2 / 

<=° -00 

x 53 Bl'm’(t') 

U',m' 

xJJ Y l>m ,(e',^)VP t (cos 1 ) dn' 


dt' (29) 
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It is clear from (29) that the vector operator (which op- 
erates on the unprimed quantities) can be moved outside 
the surface integral. Doing so, and using the following 
addition theorem of spherical harmonics [e.g., Jackson, 
1975], 


p tico *i) = YUO, t/>) (30) 

m=— t 

in (29) yields 

t,m 

X f B/(f-to}E W(0 

-io 

xJJ Y; m (0’, /)y w («', *')<«' U' (3i) 

n ' 

Finally, using the orthogonality property (5) in (31) gen- 
erates the expression 




/ - m V oo J 


xVYtmiW). 


(32) 


An analogous, though simplified, procedure applied to 
a Green’s function of the form (25) yields the following 
expression for the response: 


m^t) = 53 1 

t,m [ 


J A,(t~t')L /m (t 

—OO 

xY em (9,rp). 


y )*'| 


(33) 


The quantity in braces in (33) represents the set of 
spherical harmonic coefficients for the scalar response 
R(6,ip,t). That is, if we write the usual spherical har- 
monic expansion as 


m v-,*) = 53 RimWuo, *), (34) 

/,m 


then clearly, 


Rim(t) = f ^ A e (t - t')L im (t') dt'. (35) 


To write the vector response (32) in an analogous, purely 
spectral form requires more work. The results of the 
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appendix can be used to generate a set of coefficients 
Dimi'm' defined such that 

(6, Dtrnt’m'Ytm' («, VO • (36) 

t',m' 

For each l , m pair, only a very small number of coefficients 
Dimi'm' are nonzero. Using (36) in (32) yields, after a 
small amount of algebra, 

R(0, rP,t)=J2 R tm Y tm (0, VO, (37) 

£,m 


where 




/'.m' 


l 

J Br(t-0W(0 dt' 




(38) 


are the required spherical harmonic coefficients of the 
components of the vector response R(6, tp, t ) . 

Equations (34), (35), (37), and (38) represent the spec- 
tral expressions sought at the outset of this section. 
Equations (34) and (35) are relevant to the calculation of 
the radial displacement response and of anomalies in the 
gravitational field. Equations (37) and (38) are appropri- 
ate for analyses of the horizontal deformations associated 
with an arbitrary surface mass load. 

In the following sections we use the general results 
derived here to consider several particular applications. 
First, the three-dimensional deformational response of a 
spherically symmetric Earth, with an elastic mantle rhe- 
ology, to an arbitrary external load. Second, we general- 
ize the elastic case to consider the effects of mantle anelas- 
ticity. In both the purely elastic and anelastic cases we 
provide a spectral formalism for computing the response 
which is generalized to incorporate frequency dependent 
tidal and surface load Love numbers. In each application 
we will outline a procedure for incorporating gravitation- 
ally self-consistent ocean redistribution contributions to 
the surface mass load and consider special cases of the 
various formalisms. Finally, we use the results derived 
in this section to derive spectral response equations as- 
sociated with the surface mass load problem on a linear 
viscoelastic planetary model. 


Response to an Arbitrary External Load of a 
Spherically Symmetric Earth With an Elastic 
Mantle Rheology 

The external load to be considered in this section will 
consist of a tidal potential W(9,ip,t) and an arbitrary 
surface mass load. The direct response of the solid Earth 
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to the tidal potential is given by (6)-(8), and the response 
to the mass load is given by (32) and (33), as well as 
(34), (35), (37), and (38). In this section, the Ae{t) and 
Bt(t) will be generated from the elastic components of 
the surface load Love numbers of (11)— (13). It will be 
instructive to consider a surface mass load separated into 
severed geophysically distinct terms. For example, we can 
write 


= L A ™(0,rP,t) + L lcB (0,rP,t) 

+ L^(^0 + £ NEO (*,tM) 

+ L msc {e,xf>,t). (39) 

L A ™(9,tp, t ) is the contribution to the surface mass load 
from atmospheric pressure [e.g., van Dam and Wahr, 1987], 
L 1CE (9, ip, t) is included to represent the load associated 
with the recent retreat of small ice sheets and glaciers 
[Meier, 1984] or any mass variation associated with very 
large ice sheets such as Greenland and Antarctica (e.g., 
Zwally, 1989]. L Mlsc (0,rp,t) represents any other sur- 
face mass load not associated with ocean mass redistri- 
butions. L EO (0,ip,t) is the contribution to the surface 
mass load from gravitationally self-consistent ocean mass 
redistributions caused by the deformation of the Earth 
due to the tidal potential and the surface mass load and 
by the direct gravitational effect of these forcings. Dahlen 
[1976] and Agnew and Farrell [1978] call this component 
the “equilibrium ocean tide.” We will use the term “equi- 
librium ocean load” instead since we are including, in 
general, surface mass load contributions which are not 
associated with the Earth’s response to the tidal poten- 
tial alone. L tiEO {9, ip, t) represents the portion of the ac- 
tual ocean mass redistribution not accounted for by the 
equilibrium ocean load. L NEO (0,ip,t ) can be estimated 
by considering the difference between L EX3 {9,ip,t) and 
ocean redistributions generated using: general analytic 
treatments [e.g., Dickman, 1989], a combination of ob- 
servational constraints and theory [Schwiderski, 1978] or 
recent satellite derived observations [Christodoulidis el al, 
1988], Whereas L EO (0,xp,t ) can be computed precisely 
for a particular Earth model (see the “pseudo-spectral” 
algorithm below), L NEO (9,rp,t) is uncertain (see discus- 
sion by Dickman [1989]), motivating the partition of the 
two ocean loading terms in (39). 

For a particular surface load component L l (9,rp,t), let 
us denote the corresponding spherical harmonic coeffi- 
cients, as in (24), by L\ m . Then, using (6) and (7), (11), 
and (12) (excluding the nonelastic components), (25), 
(26), (32), (33), and (36) yields the spherical harmonic 
expansions: 


U{9,1>,t) = Y, U *mYlm{9,xP) (40) 

£,m 
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and 


(41) 


V(0,xl>,t) = '£v em (t)Y em (8,tl>), 

t,m 

where the spherical harmonic coefficients are given by 


two - 


mhr* EL ‘"'M < 42 > 


V tm {t) = 


E 

t',m‘ L 


iT,E Wif'm' (0 i 47Ta 3 iL.E r /a.s 
V — g + IWTTym: 1 * 


9 T ( 2 £' + l)M e 
xDl'm'imt 


(43) 


and 




= L?™(t)+L l £?(t) 
+ LfZ(t) + L™°(t) 
+ L™ sc (t). 


(44) 


U ( 0 , ip, t) and V ( 0 , ip, t) represent the radial and tangen- 
tial deformations due to the tidal potential W(8,ip, t) and 
the total surface mass load L(6, ip, t) acting on a spheri- 
cally symmetric planet with a purely elastic mantle. 

The various contributions to Limit) of (44) can be gen- 
erated from a spherical harmonic decomposition of the as- 
sociated space dependent fields of (39). These fields are, 
in large part, constrained by observation. For example, 
L^m U i®> Vs t) can be easily generated from meteorological 
maps of surface pressure variations. As discussed above, 
the are a notable exception to this requirement since 
they can be computed by solving an equation, known as 
the sea level equation, based on the gravitational poten- 
tial perturbation associated with the external loading of 
the planetary model. In a later section, we provide a 
straightforward pseudo-spectral algorithm for generating 
solutions of the sea level equation on a spherically sym- 
metric Earth with an elastic mantle rheology to extremely 
high spatial resolution. The algorithm outputs the coef- 
ficients Lf£ and it is thus particularly well suited for 
application within the spectral formalism for computing 
the response of the Earth described above. 

Frequency Dependent Love Numbers 


The expressions (l)-(3), for the tidal Love numbers 
and the elastic components of (11)— (13), for the surface 
load Love numbers, assume that these numbers are con- 
stant, independent of (temporal) frequency. In fact, it 
is well known that these numbers are frequency depen- 
dent, in some frequency ranges quite strongly so, and in 
this section we will generalize the results derived above 
to incorporate this particularly important aspect of the 
response of the Earth. 
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Wahr [1981a, b] has shown that the frequency depen- 
dence of the Love numbers can be expressed in terms 
of resonance expansions. Each of the Love numbers has 
a nonresonant component which is weakly dependent on 
frequency, as a consequence of the effect, on the response, 
of inertial and coriolis forces, and of the ellipticity of the 
planet. A much stronger frequency dependence is evident 
in the resonance terms of the Love number expansions 
[Wahr, 1981b] which arise from the existence of certain 
eigenfunctions in the Earth’s free response which include 
a deformational component; in particular, the Chandler 
wobble (CW), and the free core nutation (FCN), eigen- 
functions. The Love numbers at frequencies near the 
eigenfrequencies of these eigenfunctions (near 14 months 
for the CW eigenfunctions and near one sidereal day for 
the FCN eigenfunction) can differ significantly from the 
value associated with the nonresonant component alone 
[Wahr, 1981b]. The near-diumal period of the FCN eigen- 
function suggests, given the significant near-diumal com- 
ponents of the tidal forcing, that predictions of the tidal 
response of the Earth (for example) should include this 
frequency dependence. (See also Schemeck [1991].) As a 
consequence, let us write 

hl E = hJ’ B ( U ) 
lJ’ E = lJ’ E (u) 

h$’ E = h E ’ E (u;) 

{ L,E = 

k E ' B = k E ' E (w), 

where u denotes the frequency. 

To incorporate (45) into the response equations de- 
rived above, we will require explicit expressions for the 
frequency dependence of the tidal and surface mass load 
forcing. We begin with the tidal potential. 

In the analysis of Cartwright and Tayler [1971] and 
Cartwright and Edden [1973] the time dependent coeffi- 
cients, W tm {t) of (4), in the spherical harmonic expan- 
sion of the tidal potential are generated by summing a set 
of pure sinusoids of the specified frequencies and phases. 
The characteristic periods of the tidal components may 
be grouped according to the absolute value of the har- 
monic order m. For example, the m = 2 (or sectorial) 
components have a nearly semidiurnal period; the m = 1 
(or tesseral) components have a nearly diurnal period; 
and the m = 0 (or zonal) components give rise to the 
long-period tides [Melchior, 1983]. (In the vast majority 
of applications, sufficient accuracy in (4) is obtained us- 
ing degrees 2 and 3 only.) Adopting this representational 


kJ> E = kJ ,E (u) 


(45) 
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form for the time dependence in (4) gives 


N tm 

W lm {t) = J2 W? m sin (u,? m t + £? m ) , (46) 

n= 1 

where w" m , £” TO , and Wp m represent the angular fre- 
quency, phase, and amplitude, respectively, of the n th 
frequency constituent (in a total of Nt m ) of the degree £ 
and order m harmonic of the tided potential. 

In an analogous fashion, the time dependence of the 
spherical harmonic components of the surface mass load 
L( m can be expressed as 

N lm 

L tm {t) = £ L n tm sin (u£ m t + C'L) , (47) 

n= 1 


where is the phase associated with the nth frequency 
constituent of the degree £ and order m harmonic surface 
% mass load. The Ne m frequency components are included 
to incorporate all of the time dependence, at degree £ and 
order m, of both the tidal potential and the surface mass 
loads. 

Using (45)-(47), together with (36), the spherical har- 
monic coefficients (42)-(43) in the spectral response ex- 
pansions (40)-(41) may be generalized to 


N em 

Utmit) = £ 

n= 1 

47TQ 3 


A^WJ^sin + Q m ) 




+ W+l)M S {uJ ^)L7m Bin (u? m t + Am) 


(48) 


and 


V tm (t) = 

{ N e' m ' 

E 

n= 1 L 


sin (u + e?w) 


+ 


(2t + l)M ^ £ ( W ^' m ')^"' m ' S * n + Z't'm’) ^ 




(49) 


As a final point in this section, we note that not all sur- 
face mass load components are conveniently expressed in 
terms of expansions of the form (47) . One example would 
be atmospheric pressure loads with near-diumal or near- 
annual periods. These periods are sufficiently close to 
the periods of the free core nutation and Chandler wob- 
ble eigenfunctions, respectively, that the Love numbers 
which characterize the response would be frequency de- 
pendent. In these cases the forcing is more conveniently 
expressed as an integral of periodic terms rather than 
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as a sum. To consider such a forcing, the summation 
over n which appears in (47) need only be replaced by an 
integral over frequency space and (48)-(49) altered ac- 
cordingly. A second example of a forcing for which the 
form (47) is not particularly appropriate would be the 
present-day secular melting of land-based ice sheets and 
glaciers \Meier, 1984]. In this case the formalism outlined 
in the preceding section (see (40)-(44)) for the response 
to a surface mass load is most appropriate, with the load 
Love numbers chosen to reflect the long-period timescale 
of the forcing. (In this region of frequency space the Love 
numbers are very weakly dependent on frequency.) 

Solving the Sea Level Equation 

Equations (48)-(49), together with (40)-(41), provide 
expressions for the three-dimensional deformation arising 
from the tidal potential and an arbitrary surface load. 
Using (10) and the elastic component of the Green’s func- 
tion (18), together with (40) and (48), we can generate an 
analogous expression for the total gravitational potential 
perturbation at the deformed solid surface. It is 


*(W,t) = £*"(*, iM), (50) 

n=l 


where 

*"(0, rf,, t) = £ [(l + kJ' E W m ) - hJ ’ E ( W ? m )) 

£,m 

xU? m sin( W ? m t + £?m) 

+ (2<Tl)M. 0 + tf'VZn) - 

x(^ RES m (■*.«+«.) 


+i?m E ° Sill (uj m t + £"?-.)) 


W*.*) (51) 


is the nth space-time dependent frequency component of 
4>(0,t/i, t). L"^f° and £"'" m represent the amplitude and 
the phase of the nth frequency component of the spheri- 
cal harmonic coefficient (of degree £ and order m) of the 
equilibrium ocean load. The L”£ ES and are analo- 
gous terms for the total of the surface mass load contribu- 
tions other than the equilibrium ocean load. (“RES” de- 
notes the residual surface mass load.) That is, the terms 
with the superscript RES represent parameters in the ex- 
pansion (47) of the mass load L(0,ip,t) — L EO (0,xp,t) 
(see (39)) and may therefore include contributions from, 
e.g., atmospheric pressure loads, ice loads, nonequilib- 
rium ocean loads, etc. In this section we assume, re- 
gardless of their origin, that these contributions may be 
prescribed at the outset, and we will be concerned with 
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a solution for the parameters describing the equilibrium 
ocean load. 

The discussion below will refer to $ n (0,xp,t), rather 
than to the total perturbation $(8,ip, t). To simplify the 
expressions used, we define the following quantities: 

= H? m sin(u,? m t + tf m ) 

L?* es (t) = L£? ES sin(w,V + £"? m ) (52) 

Lg°(t) = L?*° S in(u,? m t + C? m ). 

Using Brun’s formula [ Heiskanen and Moritz, 1967], we 
find that the perturbed position of the equipotential sur- 
face (originally located at the undeformed solid surface) 
with respect to the deformed solid surface, given the 
gravitational potential perturbation (51), is 4> n ($, xp, t)/g. 
This is not, however, the location of the perturbed sea 
level, or geoidal, surface with respect to the deformed 
solid surface since that surface is not constrained to re- 
main on the same equipotential. The location of the per- 
turbed geoidal surface with respect to the deformed solid 
surface can be written as <f> n ($, xp, t)/g-\- A$ n (t)/g, where 
A $>"(<) is a time dependent quantity whose value is con- 
stant over the surface of the Earth. The perturbation in 
the ocean thickness, or bathymetry, is simply the projec- 
tion of this field onto the region of the Earth covered by 
oceans. 

Let us define the ocean function, C(6,xp), as a field 
equal to unity over the oceans, and zero elsewhere [Munk 
and MacDonald, 1960], and C/ m as its set of spherical 
harmonic coefficients. Let us also define the equilibrium 
ocean load thickness by S EO (0, xp, t) and its associated 
coefficients by Sf£(t). (lF°{0 y xl),t) = p w S EO (0,xp,t), 
L /m W = Px«S^{t), and L”£° (t) = where 

p w is the density of sea water.) Then, using the argu- 
ments above, the equation which governs Sj*|f°(0,V»,t) 
is 


£s£f°(t)y /m (0,vo « r,Y rt (o,xp) 


£(% 

/,m ' 




AS> n (t) 

9 


SioSmOj Yi, 


.(*,*) 


, (53) 


where the are the coefficients of the spherical har- 
monic expansion of the total gravitational potential per- 
turbation $ n (6,xp,t). These coefficients are provided in 
(51). They are 


= (l + kJ’ E W, n) - hJ’ E (w? m )) WS*(t) 
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Equation (53) is the sea level equation governing the 
equilibrium ocean load height on a spherically symmetric 
Earth with an elastic mantle rheology subject to external 
loading. The equation derives from an integral equation 
for the sea level response, manifested in the appearance of 
the coefficients °(t) on both sides of the equation via 
(54). Physically, this equation reflects the fact that the 
equilibrium ocean load perturbs the gravitational field of 
the planet and in turn is governed by this perturbed field, 
since the geoid is defined to be an equipotential surface. 

An expression for the quantity A $ n (f) in (53) may be 
obtained by imposing a conservation of mass constraint 
on the equilibrium ocean load. Integrating both sides of 
(53) over the surface area fi yields (using the orthogo- 
nality property (5), and the fact that Yoo = 1 in our 
normalization), 


A *”(Q S£j EO (t) 1 1 

g Coo Coo 


*[[ dsi. (55) 

n r,« t,m 


One may directly obtain an expression for S^q EO (£) by 
considering all of the surface mass load contributions in 
(44) which are reflected in the mass balance of the global 
oceans. For example, if an application were being con- 
sidered in which only L ICE (0, ip, t) affected this mass bal- 
ance, then imposing a conservation of mass yields 

SS ‘6 EO W = ~L$ CE (t). - (56) 

Pw 

Dahlen [1976] has provided a purely spectral algorithm 
for solving a somewhat less general, but essentially equiv- 
alent equation of the form of (53). The algorithm involves 
multiplying (53) by a basis set Yp q {6,ip), integrating over 
fl, and making use of the orthogonality property (5). The 
procedure requires the determination of coupling coeffi- 
cients between sets of surface spherical harmonics, and 
this severely limits the maximum degree and order (that 
is, truncation level) which can be feasibly considered. 
( Mitrovica and Peltier [1991a] have argued that a trun- 
cation at degree and order near 32 is now the maximum 
currently feasible.) A hybrid spectral technique is neces- 
sary to compute the equilibrium ocean load at the higher 
resolution required, for example, near continental mar- 
gins. We outline below a straightforward pseudo-spectral 
approach which can be applied with high truncation lev- 
els. (We have found that consideration of degrees and 
orders up to 1024 is certainly feasible.) 

Let us represent an approximate solution to the equi- 
librium ocean load height by {S^f°(t)}'. We can use 
(53) iteratively to improve the estimate of (t) by 
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substituting {S7rrf° (t)}* into the right-hand side of the 
equation. This substitution yields 


= cm) 

£,m 


X 



,(57) 


where the superscript t+1 denotes the improved estimate, 
and, using (54), 


*5f «) - [i + 

+ VF?i j k [‘ + ^ E (“S,)-^ E W] 

X (l£j res (0 + (*)}’) • (58) 

Let us define 

V' n M,t) = 5^^(t)y P ,(0,t&), (59) 


and, furthermore, the projection of this field onto the 
ocean regions as 

(*,<M) = ^•”(</,V’,f)C(tf,V'), (60) 

which has a spherical harmonic expansion given by 

*& E (*,tM) = ^ {^c E (0} fm ^m(«,V’). (61) 


Using (60) and (61) in (57) yields, for each spherical har- 
monic coefficient, 


(62) 

Isolating the degree and order zero term in (62) gives 


A$ i,n (t) 

9 



(63) 


Equation (63) can also be obtained by applying the im- 
provement procedure to (55). The term Sfy EO (t) does 
not need to be iteratively improved upon since it can be 
determined a priori from the spectral coefficients of the 
appropriate surface mass load contributions (as in (56), 
for example). 

Equations (58)-(63) define the pseudo-spectral algo- 
rithm for computing the equilibrium ocean load height 
on an externally loaded spherically symmetric Earth with 
an elastic mantle rheology. The procedure begins with an 
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approximation to the spherical harmonic coefficients of 
the equilibrium ocean load, {S£^ EO (t)}\ Using (58), the 
response coefficients $>)>£ (£) can then be determined, and 
the field $ t,n (6,rp,t) of (59) can be synthesized at spe- 
cific points on a latitude-longitude grid over the Earth’s 
surface, with the grid spacing dependent on the spheri- 
cal harmonic truncation level. These grid values are then 
multiplied by the associated value of the ocean function 
(that is, either 0 or 1; (60)), and the resultant grid de- 
composed to yield the {$oceW}<”>- These coefficients 
are used in (63) to determine A $*’"(£) and then, to- 
gether with this value, in (62) to generate the next, or im- 
proved, estimate {•Sj^f°(£)}* +1 . The procedure is termed 
pseudo-spectral because operations are performed in both 
the spectral and space domain, which is in contrast to the 
Galerkin approach of Dahlen [1976], 

The starting value {S" r ^ EO (t)}‘ =1 will depend on the 
particular application of the procedure (see below) . As a 
convergence parameter one may compute 




(64) 


We have found, in applications involving a variety of ex- 
ternal forcings (including tidal forcings, atmospheric pres- 
sure loads and ice loads), that three iterations are suffi- 
cient to ensure that £* +1 < 10 -4 . 

The formalism outlined above is related to the pseudo- 
spectral approach of Mitrovica and Peltier [1991a] for com- 
puting the gravitationally self-consistent redistribution of 
ocean mass arising from the melting of late Pleistocene 
glaciers on a linear viscoelastic Earth model. The proce- 
dure described here, however, considers the elastic case 
and generalizes the formalism to include a gravitational 
potential forcing as well as a variety of surface mass loads, 
which may or may not affect the mass balance of the 
global oceans. The present formalism also incorporates 
frequency dependent Love numbers. 

As derived above, the pseudo-spectral formalism is ap- 
plied to each frequency component of the external forcing. 
In the event that a frequency range is being considered 
within which the Love numbers do not change apprecia- 
bly, then the formalism can be applied once to generate 
the response over that frequency range. In this case, the 
governing equations can be derived by suppressing the 
frequency dependence of the Love numbers in the equa- 
tions above, and by considering the “total” (rather than 
the frequency dependent) harmonics in the specified full 
frequency range. That is, all terms of the form 
can be replaced by £e m [t), where it is understood that 
the term is equal to the sum of the £? m (t) over 

the appropriate range of n. The same algorithm is appro- 
priate when considering the equilibrium ocean response 
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arising from a surface mass and having a time depen- 
dence which is secular, rather than periodic, in nature. 
In this case the Love numbers are chosen as those appro- 
priate to very low frequencies, and the response and load 
harmonics in the above equations refer to the total fields 
rather than any specific frequency constituent (or range 
of constituents). 

Below, we briefly consider two special cases of the gen- 
eral formalism outlined above. 

Self-Consistent Equilibrium Ocean Tides 

Here we consider the equilibrium ocean load height 
arising solely from the application of a zonal (that is 
m = 0) tidal potential forcing of degree 2. Such a forc- 
ing gives rise to the zonal, or long-time scale, body tides 
[Melchior, 1983], and the general equation (4) for the tidal 
potential can be replaced by 

W(0,iJ;,t) = W m (t)Y2o(0,1>). (65) 

Using (46) and (52), we have 

Nio 

WaoW = ^"o g i n (wzo* + £20) 

n=l 

N20 

= ( 66 ) 

n=l 

In this case, the response coefficients $>" m (£) in the sea 
level (53), given by (54), are replaced by 

*X»W = (l + w? m (t)6„6 m0 

+ x& r i % 0 + - h$’ E wj) 

(67) 

Equation (58) is changed accordingly. Furthermore, since 
the mass of the oceans remains constant in this case, 
5- eo (£) = 0 in (55), and (63) becomes 

A$ i,w (t) l fSocEffl}^ 

9 9 Coo 

The pseudo-spectral formalism proceeds as above (equa- 
tions (58)-(62) and (68)) and can be initiated using 

{s£f°M}‘ =1 = o- 

Self-consistent equilibrium ocean tides induced by a 
zonal tided potential have edso been considered by Mer- 
riam [1973, 1980] and Agnew and Farrell [1978] in their 
analysis of the effect of such tides on the length of day. 
We have compared ocean tide results generated using the 
pseudo-spectral approach with results presented for a spe- 
cific W 20 amplitude by Agnew and Farrell [1978, Figure 1] 
and have found qualitative agreement. 
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Ocean Redistributions Due to Present-Day Ice Melting 

Here we consider the equilibrium ocean mass redistri- 
bution arising from the melting of present-day land-based 
ice sheets and glaciers \Meier, 1984]. In this case, (54) and 
(58) reduce to 

= (2e+\?M e (* + k t' E ~ h *’ E ) 

( 6Q ) 

and 

*!"•<*> " WTT?m; (‘ + k ‘ E - k ‘‘ B ) 

xUCW + p44 0 W}‘), (70) 

respectively, where pi is the density of ice and where 
p//]^ E (t) = L\^(t). represent the time depen- 

dent spherical harmonic coefficients of global ice sheet 
thickness. Equation (56) thus becomes 

S 0 E o°(0 = /£ E (0, (71) 

Pw 

which can be applied to (63). 

A good first guess to the ocean redistribution is the 
so-called eustatic redistribution. That is, 

{^(t)} <=1 = (72) 

Pw Coo 

Equation (72) represents the case of meltwater being re- 
distributed in the oceans in a manner independent of ge- 
ographic location. 

As discussed above, in this case, where the time de- 
pendence of the surface mass load is secular, rather than 
periodic, the pseudo-spectral formalism is applied to com- 
pute the total response rather than a specific frequency of 
that response. (The superscript n has been dropped from 
the equations.) The Love numbers in (69) and (70) refer 
to those appropriate for a long-period surface loading. 

Incorporating an Anelastic Component of 
Planetary Response 

The Love numbers described above were associated 
with an Earth model having a purely elastic mantle rheol- 
ogy. The effects of mantle anelasticity on the Earth’s re- 
sponse to relatively short timescale external forcings can 
be incorporated into these Love numbers using a straight- 
forward numerical perturbation technique described in 
detail by Wahr and Bergen [1986]. The Wahr and Bergen 
[1986] approach yields complex Love numbers which re- 
place the real coefficients (l)-(3) and (the elastic com- 
ponents of ) (H)— (13), reflecting the fact that anelas- 
ticity perturbs the amplitude of the Love numbers and 
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introduces an out-of-phase component (or “lag” ) in the 
response. 

Let us denote a complex Love number of arbitrary type 
(i.e., tidal or load, and h, k, or l) and degree as S, where 

S=S / + tS°, (73) 

and I and O denote the in-phase and out-of-phase com- 
ponents, respectively, of the Love number. Let us also 
assume an external forcing (that is, an external potential 
or a surface mass load) of the form Fcos(u/t -f 4>), and let 
us construct from this a generalized function of the form 

F = Fe'ivt+t)' ( 74 ) 

The response of the planet to the forcing is given by 
H = Re{5 • F} 

= S 1 F cos(ut + <f>) — S°Fsin(iJt + <f>), (75) 

where Re denotes the real part of the complex quantity. 
The first term on the right-hand side of (75) is the in- 
phase component of the response, whereas the second 
term is the out-of-phase component of the response which 
arises solely in consequence of anelasticity. «S° can be 
computed for any particular model for mantle anelastic- 
ity using the technique of Wahr and Bergen [1986]. 

The form (75) is introduced to emphasize that the out- 
of-phase component of the generalized response, the sine 
term of (75), can be computed using the spectral formal- 
ism outlined above for the elastic response. One needs 
only to replace the external forcing terms by their out- 
of- phase analogues and the Love numbers by their imag- 
inary components. The spectral calculations of three- 
dimensional crustal deformations, equilibrium sea level 
variations, etc., then proceeds in entirely the same fash- 
ion. 

For external forcings on longer timescales (for exam- 
ple, the surface mass load associated with the melting of 
the late Pleistocene ice sheets), the nonelastic behavior 
of the planet is evident in a viscoelastic response of the 
planet [e.g., Cathles, 1971; Peltier, 1974]. In this case it is 
appropriate to consider a viscoelastic formalism for the 
Earth’s response, as we do below. 

Response to an Arbitrary Surface Load of a 
Spherically Symmetric Earth with a 
Viscoelastic Mantle Rheology 

In this section we consider three-dimensional deform ar- 
tions associated with the surface loading of a spherically 
symmetric Earth having a linear viscoelastic mantle rhe- 
ology. In this case we need only consider the response (32) 
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and (33) with the appropriate Green’s functions. Using 
(12) and (15) in (32), with Bi defined as in (26), yields 


U(W,t) = 


E 

/,m 


K(t) 


(2€ + 1 )M e I + 


fc=l 


xY em (6,rP) 


(76) 


and 

«'(»,*<) - 


E 

/,m 


K(t) 


( 2 riiw; Llm{t)l *' E + £ r2 *’^ 


*= i 


x^y/ m (0,v>), 


(77) 


where 


/£»(*) = t Limit') exp [~4(t - t')] dt>. (78) 

«/ — oo 


The first terms in large parentheses on the right-hand 
sides of (76) and (77) represent the purely elastic compo- 
nents of the response, whereas the second terms represent 
the nonelastic response. 

Predictions of present-day three-dimensional deforma- 
tions associated with the melting of the massive late Pleis- 
tocene ice sheets represent an obvious application of the 
above theory. It has been common to model the time 
dependence of these ice sheets as a series of step-load in- 
crements [Farrell and Clark, 1976; Wu and Peltier, 1983]. 
This representation incurs no loss of generality since the 
time steps may be made arbitrarily small. If we apply the 
same discretization to our general surface load L(0,ip,t), 
we can write 

N 

L{0,ip,t) = Y,6L n (0,rP)H(t-t n ) 

n= 1 
N 

= ( 79 ) 

n=l l,m 

where H represents the Heaviside step function and N 
is the total number of step load increments. Comparing 
(79) with (24) gives 


N 

W*) = ][>£*»//(*-*„). ( 80 ) 

n=l 

The form of (80) allows for an analytic evaluation of the 
integral in (78). Applying this form to (76) and (77) 
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yields, using (36), response equations of the form (40)- 
(41) with spherical harmonic coefficients given by 


U tm (t) = 


47ra 3 


(2 If + 1 )M, 




Ltt \ 

+ SS^ L "m !1 4-( 1 - eX P[-4( t - t n)3) (81) 

*=ln=l S * / 


and 


*<OW< 0 r0 LS , . 1 -\'| 

+ £ J2 6L t'm'-$r- (l - exp [-«£(t - «„)]) ) [ 

*=1 n=l Sk / ) 


^ Dfm'tm • 


(82) 


The summations over n include all the load increments 
which have a time onset prior to the observation time t. 
That is, N(t) represents the number of load increments, 
in the total of N, which satisfy t„ < t. 

As above the surface mass load appearing in (81)-(82) 
can include contributions from a number of different geo- 
physical loading phenomena acting on a broad range of 
timescales (see (39)). Following (39) and (44) we can 
separate the load into a contribution from a gravitation- 
ally self-consistent equilibrium ocean load and a residual 
( “RES”) surface load which includes all other surface load 
contributions. That is, we can write 

Limit) = PvSffiit) + lgf(<) (83) 

and 

6L n tm = p w 6S^° + 6L^ RES , (84) 

where the S^(£) are the spherical harmonic coefficients 

_ DA 

of the equilibrium ocean load height and the tire 

the associated coefficients for the nth incremental change, 
as expressed in (80). 

Mitrovica and Peltier [1991a] have outlined both a spec- 
tral and a pseudo-spectral formalism for computing grav- 
itationally self-consistent equilibrium ocean load height 
increments associated with the surface loading of a spher- 
ically symmetric Earth with a viscoelastic mantle rheol- 
ogy. Both algorithms output the spherical harmonic coef- 
ficients of the equilibrium ocean load height (or sea level) 
increments, f° and are thus ideally suited (as a pre- 
liminary step) to the spectral formalism outlined here for 
computing the three-dimensional deformations associated 
with an arbitrary surface load. 
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Present-Day Deformations Associated With the Melting 
of the Late Pleistocene Ice Sheet 

Predictions of present-day deformations associated with 
the melting of the late Pleistocene ice sheets can be based 
upon applications of (81)-(82) using quite specific forms 
for (83) and (84). These forms are 

Limit) = p w S™(t) + P/ /J£ E (t) (85) 

and 

6L?m = P W 6S?*° + PI SI^ CE , (86) 

where /]^ E (t) are the spherical harmonic coefficients of 

ipp 

the late Pleistocene ice sheet thickness, and 6iy^ are 
the associated coefficients for the n th incremental change 
(as expressed in (80)) in ice thickness. The coefficients 
for Jj£ E (t) ^d 6 ^E can be generated from published 
late Pleistocene ice thickness reconstructions [e.g., Tush- 
ingham and Peltier, 1991], while the 5^(t) and 
can be determined using the pseudo-spectral algorithm 
developed by Mitrovica and Peltier [1991a]. 

Perturbations in the Center of Mass of the 
Deformed Planet: A Spectral Approach 

Let us use the term “planet” to include the Garth 
model, whose undeformed state is used in the calcular 
tion of the Love numbers, but not in the calculation of 
the surface mass load. In this section we derive, using 
a spectral formalism, an expression for the perturbation 
in the center of mass of the planet, relative to the cen- 
ter of mass of the undeformed state, which arises as a 
consequence of a surface mass loading. 

Cathles [1971] has recognized that the center of mass 
of the planet can be perturbed in the surface mass load 
problem, since the load may have a nonzero spherical har- 
monic degree one forcing. We will presume in this section 
that no other external forces (e.g., tided) are acting on the 
system. Accordingly, perturbations in the center of mass 
of the planet and surface load are related by 

Arf“ _ m 

where r° M denotes vectors measured relative to the un- 
perturbed center of mass of the system, e and SL denote 
the planet and surface load, respectively, Ms l is the mass 
of the surface load, and A denotes a perturbation. Since 
the surface load conserves mass, we can write 

Ar s c L M = JJJv A p(r, 6, ip) r dv, (88) 

where Ag is the position dependent density perturbation 
in the surface mass load and the integration is performed 
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over the volume of the load. The confinement of the sur- 
face mass load to a thin surface layer allows us to write 


A e(r, 6, VO = Acx(0, VO 6(r - o), (89) 


where 6 is the delta function and A<r is the perturbar 
tion in the surface mass load, in units of mass per area 
in accord with our definition of the surface load used in 
previous sections. Using (89) in (88) yields 


Ar<£ M = A<r(0,VO f sin 0 d6 dip, 

n 


(90) 


where f is the unit radial vector. If we denote the spheri- 
cal harmonic coefficients of the surface mass load pertur- 
bation A a{0,ip) by Aa> m , then the surface integration in 
(90) can be shown to be 


JJ Aa(0,rp) f sinfl d9 dip = Acri m , (91) 

n m=_1 


where the e m are the polarization unit vectors defined by 
Edmonds [1960] and given in the appendix, and Acri m are 
the coefficients in the spherical harmonic expansion of the 
surface mass load. Combining (87), (90), and (91) yields 


Ar e CM = - 


47ra 3 
M e y / 3 



(92) 


If, instead of deformations, we are interested in deforma- 
tion rates, then we can write 


Ar, 


CM 


47Tfl 3 

M c yJ 3 


1 

^ ^ Cm Adj m , 

m =— 1 


(93) 


where the A<7i m are the degree one coefficients in the 
spherical harmonic expansion of the time rate of change 
of the surface mass load. The vector components of the 
center of mass perturbations given by (92) and (93) are, 
using the polarization unit vectors listed in the appendix, 
referenced to the geocentric Cartesian coordinate system. 

The Love numbers for an Earth model are computed 
so that the Earth response generated from them is re- 
ferred to an origin at the center of mass of the undeformed 
planet [Farrell, 1972]. Thus positions computed losing the 
spectral response equations derived above will be refer- 
enced to this origin. Equations (92) and (93) provide 
the perturbation (or the rigid translation) of the center 
of mass of the deformed planet with respect to its un- 
deformed state, and they may therefore be used in an a 
posteriori calculation to reference any computed changes 
in position to an origin at the center of mass of the de- 
formed planet. In other words, (92) and (93) may be 
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used to remove the component of pure translation asso- 
ciated with position changes computed in reference to an 
origin at the center of mass of the undeformed planet. 
The positions thus corrected are more accurately termed 
“deformations.” 

Summary 

We have outlined a complete spectral formalism for 
computing three-dimensional deformations of a spheri- 
cally symmetric Earth due to an arbitrary surface mass 
loading. The appropriate spectral response equations 
have been derived for the case of Earth models with elas- 
tic, anelastic, and viscoelastic mantle rheologies. Method- 
ologies required for computing the Love numbers appro- 
priate for such Earth models are found in the literature; 
the required references are Farrell [1972], Wahr and Bergen 
[1986], and Peltier [1974, 1985], respectively. In the case of 
an elastic or anelastic mantle rheology, we have incorpo- 
rated frequency dependent Love numbers into the formu- 
lation and have derived a pseudo-spectral algorithm for 
determining the ocean mass redistribution arising from 
the gravitational perturbation associated with the exter- 
nal loading. As a consequence of their definition, response 
predictions based on tidal and surface load Love numbers 
are generally referenced to an origin located at the cen- 
ter of mass of the undeformed planet. Coordinate po- 
sition changes may thus incorporate a pure translation 
associated with the perturbation of the center of mass of 
the deformed planet from its location in the unperturbed 
state. Accordingly, we have derived a spectral formula 
for translating the origin to the center of mass of the de- 
formed planet, and the position changes thus “corrected” 
are more accurately interpreted as “deformations.” 

In order to test numerical codes established to apply 
the theory described herein, we have compared our pre- 
dictions of three-dimensional crustal deformations with 
published results based on the Green’s function approach 
described in the introduction. Those tests indicate that 
the spectral formalism can be easily applied to spherical 
harmonic truncation levels in excess of degree 2048 using 
existing computer systems. 

In a companion article [Mitrovica et al„ this issue] we 
generate predictions, up to degree 512, of present-day 
three-dimensional deformations associated with the glac- 
ial isostatic adjustment process. The deformation pre- 
dictions are based on the spectral response equations de- 
rived above, and they incorporate both gravitationally 
self-consistent ocean redistributions (using the pseudo- 
spectral algorithm of Mitrovica and Peltier [1991a]) and 
realistic models of the space-time deglaciation history of 
the late Pleistocene ice sheets. 
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Appendix: An Analytic Form for VY* m (i9, ip) 


Forte and Peltier [1994] have derived the following iden- 
tity for VYt m (0,ip) where the two-dimensional gradient 
operator V is defined as in (9) and the surface spherical 
harmonics are normalized as in (5): 


vy /m (0,VO = 


e_i 

y/2 


1W, («,*)< ( (<+ (S+lii2/ + " , 3) +1> ) 
+ 14-.W*. *)«+« ( ( fo+I)$-~i)‘ ) ) 
*«.-<*•*>* 

(/ + i) (§rrpni) 

+ %[ Y ‘ ♦■—■ <«•*>' 


(Al) 


where e_j, &o, and ei are the polarization unit vectors 
defined by Edmonds [1960] as 

e-i = -j=(x-iy) (A2) 


and 


II 

o 

<v 

(A3) 

ei = —j=(x + ty). 

(A4) 


The vectors x, y and z represent the basis unit vectors in 
the geocentric Cartesian coordinate system. The vector 
YYi m (0, VO can, of course, also be written in the form 


Vr* m (0, VO = x X + yY + zZ, (A5) 


where the components X, Y, and Z can be determined 
by substitution of (A2)-(A4) into (Al). It will be useful 
to generalize (Al) by writing it as 

OO t 

Vy, ro (0,VO=£ Yj D tmt >m'Yr m >(0,rl>). (A6) 

t'=0m’=-t' 


In order to express VY* m (0, VO in terms of a particular 
basis set, it is necessary only to determine the vector 
coefficients Dtmi'm', which, using (A1)-(A4), amounts 
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to determining expressions for the basis set of interest in 
terms of the (e_i,eo,ei) or the (x, y, z). It is clear from 
(Al) that, for a given £, m pair, only a very small number 
of Derntm' are nonzero. 

The surface spherical harmonics Yt m (0, ip) have the fol- 
lowing property: 


Ip) = (-iry/ m (0, IP), (A7) 

where the asterisk denotes the complex conjugate. One 
can show, using (A7) in (A1)-(A4), that the vector com- 
ponents of VY tm (6, ip) have the analogous property: 

VY e ,- m (0,ip) = (-i) m [vy fm (tf,v>)]*. (AS) 

Let us consider a vector defined such that 

oo £ 

W,tP)=Y l ^2 T lm VY tm (e,rP). (A9) 

e=0 m=-/ 

If the components of T(0, ip) are real, then one can show, 
using (A8), that 


Tl _ m = (— l) m T/ m , (A10) 

and that the summation (A9) may be written as 

oo / 

£=Om=l 

oo 

+52TtoVY'o(0,<P), (All) 

/=o 

where Re denotes the real part of the complex quantity. 
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Abstract. Using a spherically symmetric, self-gravitating, linear viscoelastic Earth 
model, we predict present-day three-dimensional surface deformation rates and 
baseline evolutions arising as a consequence of the late Pleistocene glacial cycles. 
In general, we use realistic models for the space-time geometry of the final late 
Pleistocene deglaciation event and incorporate a gravitationally self-consistent ocean 
meltwater redistribution. The predictions of horizontal velocity presented herein 
differ significantly, in both their amplitude and their spatial variation, from those 
presented in an earlier analysis of others which adopted simplified models of both 
the late Pleistocene ice history and the Earth rheology. An important characteristic 
of our predicted velocity fields is that the melting of the Laurentide ice sheet over 
Canada is capable of contributing appreciably to the adjustment in Europe. The 
sensitivity of the predictions to variations in mantle rheology is investigated by 
considering a number of different Earth models, and by computing appropriate 
Frechet kernels. These calculations suggest that the sensitivity of the deformations 
to the Earth’s rheology is significant and strongly dependent on the location of 
the site relative to the ancient ice sheet. The effects on the predictions of three- 
dimensional deformation rates of altering the ice history or adopting approximate 
models for the ocean meltwater redistribution have also been considered and found 
to be important (the former especially so). Finally, for a suite of Earth models we 
provide predictions of the velocity of a number of baselines in North America and 
Europe. We find that, in general, both radial and tangential motions contribute 
significantly to baseline length changes, and that these contributions are a strong 
function of the Earth model. We have, furthermore, found a set of Earth models 
which, together with the ICE-3G deglaciation chronology, produce predictions of 
baseline length changes that are consistent with very long baseline interferometry 
measurements of baselines within Europe. 


Introduction 

The last late Pleistocene deglaciation event of the cur- 
rent ice age was sufficiently massive and recent that the 
Earth remains in a state of appreciable isostatic dis- 
equilibrium. This disequilibrium is manifest in a va- 
riety of geophysical observables, but none more direct 
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than the associated three-dimensional crustal deform a- 
tions and velocities. In this regard, radial deformations 
have received the most consideration, mainly because of 
their important contribution to apparent sea level variar 
tions in both the near field and the far field of the late 
Pleistocene ice sheets [Cal hies, 1975; Peltier, 1986; Mitro- 
vica and Peltier, 1991]. In spite of this interest, however, a 
number of fundamental questions regarding radial crustal 
deformations remain unanswered. In particular, changes 
in the depth-dependent sensitivity of the radial velocity 
field to variations in the mantle viscosity profile, as one 
proceeds from sites in the near field of the late Pleistocene 
ice sheets, through the peripheral bulge, and on to the far 
field, has not been investigated. A small set of forward 
calculations, based on spherically symmetric Earth mod- 
els having isoviscous upper and lower mantle regions [e.g., 
Cathles, 1975; Peltier, 1986], suggests that these changes 
must be significant. 

The study of James and Morgan (1990] represents the 
first to consider horizontal motions associated with the 
glacial isostatic adjustment phenomenon. James and Mor- 
gan [1990] considered, in particular, present-day horizon- 
tal deformations and velocities arising from the melting 
of a Laurentide scale ice sheet; their analysis included 
a number of significant approximations. The Laurentide 
ice sheet was modeled as a single ice disk having a fixed 
basal radius and a mass which varied, in time, with a 
simple saw-tooth pattern. The Earth model adopted in 
the calculations was characterized by a five-layer, incom- 
pressible rheology. The tangential velocity profiles com- 
puted by James and Morgan [1990] were characterized by 
a broad region, extending from the center of the ice sheet 
to well into the far field, in which tangential motions were 
directed inward (i.e., toward the disk center). A recent 
investigation by James and Lambert (1993] suggests that 
the James and Morgan [1990] predictions do not closely 
approximate those obtained using more realistic models 
of the late Pleistocene surface mass load and the Earth’s 
rheology. 

Current and future space geodetic measurements of 
relative three-dimensional crustal deformation rates tire, 
and will continue to be, of sufficient accuracy to require a 
comprehensive examination of present-day three-dimens- 
ional deformation rates associated with the glacial iso- 
static adjustment phenomenon. This paper provides such 
an analysis. As a preliminary, we examine predictions 
generated using simplified disk models of the late Pleis- 
tocene ice sheets, our goal being to isolate some of the 
basic characteristics and sensitivities of the postglacial 
velocity field. Our main predictions will, however, be 
based on a surface mass load composed of more realistic 
ice and ocean load components. 

The sensitivity of the predictions to the radial profile 
of mantle viscosity will be examined by using a suite of 
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forward calculations and by computing Frechet kernels. 
In particular, we compute the kernels for predictions of 
present-day radial and tangential velocities at a number 
of sites in North America and Fennoscandia which are 
chosen to sample the near, intermediate, and far field 
glacial isostatic adjustment process. In general, the calcu- 
lations will be performed by adopting the ICE-3G model 
of the late Pleistocene deglaciation event; however, we 
consider the sensitivity of the forward predictions to vari- 
ations in the ice model. 

Finally, we also present predictions of the evolution of 
a number of North American and European baselines in- 
cluded in a recent analysis of Mark III geodetic very long 
baseline interferometry (VLBI) data [/tyon et al., 1993]. 
We have several goals: (1) to obtain values for the po- 
tential range of signals in the baseline evolution from the 
glacial isostatic adjustment process, (2) to consider the 
sensitivity of the predictions to variations in the mantle 
viscosity profile, (3) to investigate the relative contribu- 
tions to the components of the evolution associated with 
radial and tangential motions of the sites defining the 
baselines, and (4) to use the predictions to investigate 
whether the VLBI-determined European baseline rates 
may be fit by deformations associated with any of a num- 
ber of “reasonable” combinations of ice history and Earth 
rheology. 

Results 

Mitrovica et al [this issue] outline a spectral formal- 
ism for computing three-dimensional deformation rates 
on a spherically symmetric, self-gravitating, Maxwell vis- 
coelastic planetary model due to an arbitrary surface 
load. We have adopted this formalism to predict glacial 
isostatic adjustment rates using a truncation at spheri- 
cal harmonic degree and order 512. The redistribution of 
mass in the global oceans, when it is incorporated in the 
loading calculations, is computed using the gravitation- 
ally self-consistent pseudo-spectral algorithm of Mitrovica 
and Peltier [1991]. In the following two sections we con- 
sider the three components of velocity in the local carte- 
sian reference frame, (local east), U$ (local south), 
and U T (local up, or radial). The three-dimensional base- 
line velocity vectors described in the final section are 
expressed in terms of components in an Earth-centered 
baseline-oriented coordinate system defined by the base- 
line component axes. The three baseline component axes 
are defined by the unit vectors l (“length”), i (“trans- 
verse"), and v (“vertical”), specified by 
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where fj and r 2 are the a priori site position vectors. The 
length, transverse, and vertical velocity components will 
often be referred to below as L, T, and V, respectively. 
Radial site motions contribute only to the L and V com- 
ponents of baseline changes, whereas tangential motions 
will contribute to all three components. 

The numerical algorithm used to determine Frechet 
kernels is described by Mitrovica and Peltier (1993a). These 
kernels relate arbitrary perturbations in the radial viscos- 
ity profile to the consequent perturbation in a particular 
observable. If we denote the observable by U, and its 
perturbation by 6U, then these kernels, which we denote 
by F(i/(r), r), are defined via the following relationship: 


l 

6U = J F(u(r), r) <5 log v(r) dr (2) 

r CMB 

where r denotes the nondimensional radius (i.e., normal- 
ized by the radius of the Earth), i/(r) is the radial vis- 
cosity profile, and tcmb is the (normalized) radius of the 
core-mantle boundary (CMB). 

Finally, all calculations will use an Earth model char- 
acterized by an elastic structure given by the seismically 
determined Preliminary Reference Earth Model (PREM) 
[Dziewonski and Anderson, 1981], an inviscid core, an elas- 
tic lithosphere (of thickness denoted by LT), and isovis- 
cous upper and lower mantle regions (with viscosities de- 
noted by i sum and vlm, respectively). 

Calculations Using Simplified Disk Load Models 

In this section we present calculations of present-day 
three-dimensional deformation rates obtained using two 
simplified models of ice loading. No ocean load will be 
included. The ice models consist of single disks having 
parabolic vertical cross section and circular horizontal 
cross section. At glacial maximum the two disks have 
basal radii of 16° and 8° and provide a rough approxi- 
mation to the maximum geographic extent of the Lau- 
rentide and Fennoscandian ice sheets, respectively. The 
time dependence of the disk mass is a simple sawtooth 
function with a 90-kyr growth phase followed immedi- 
ately by a 10-kyr deglaciation phase. Within each cycle 
the radius and height of the ice disk are varied in order 
to maintain plastic equilibrium ( Paterson , 1981). The final 
deglaciation event for the ice disk models will be assumed 
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to have ended 6 kyr B.P. The Fennoscandian ice sheet dis- 
appeared earlier than this date; however, we will be most 
concerned in this section with isolating the influence of 
the horizontal scale of the ice sheet on the numerical pre- 
dictions. The maximum height of the ice disk is chosen to 
yield a peak present-day uplift rate of 10 mm yr -1 when 
we use the “test” Earth model described below. 

In Figure 1 we show results obtained using the two 
ice disk models and a suite of Earth models. In this and 
subsequent sections we will refer to a “test” Earth model. 
This model is defined by LT = 120 km, uum = 10 21 Pa s, 
and ulm = 2 x 10 21 Pa s and has been proposed by 
Tushingham and Peltier [1992] on the basis of an analysis of 
a globally distributed database of relative sea level (RSL) 
curves. The Earth models considered on Figure 1 include 
the test model and others identical to it except that one 
of either LT, vum, ° r t>LM has been varied. 

To begin, let us consider the results for the Laurentide 
scale ice sheet. As is well known, the radial velocity field 
is characterized by a region of uplift, or rebound, over the 
bulk of the previously glaciated region, encircled by a pe- 
ripheral region undergoing subsidence. The pattern of ra- 
dial motions associated with the deglaciation is relatively 
insensitive to variations in lithospheric thickness. Within 
the central region the uplift rates diminish as either uim 
or uum is reduced: a reduction in the viscosity (over the 
range considered in Figure 1) leads to Earth models with 
shorter characteristic relaxation times; hence these mod- 
els will have adjusted, by the present, to a state closer 
to their final equilibrium position. As vlm is increased, 
the position of maximum subsidence within the periph- 
eral bulge migrates away from the previously glaciated 
region. Within the bulge region the most pronounced ef- 
fect of increasing v\jm i8 to increase the subsidence rate. 

Next, consider the tangential velocity profile associated 
with the test model. In the region extending from the cen- 
ter of the now vanished disk to the location of maximum 
present-day peripheral subsidence rate, the tangential ve- 
locities are directed away from the disk center. Beyond 
this central region the “test” profile indicates a very ex- 
tensive region in which the tangential velocity is directed 
toward the previously glaciated region. Indeed, one must 
move approximately 70° from the disk load center before 
the amplitude of this negative velocity component is re- 
duced by 50% of its peak value of —0.3 mm yr - 1 . This re- 
sult suggests that the melting of the Laurentide ice sheet 
has the potential to produce appreciable present-day tan- 
gential velocities in regions as distant as Fennoscandia 
(displaced only about 45° from Hudson Bay) and, indeed, 
well beyond. 

The tangential velocity field is sensitive to variations 
in the lithospheric thickness. The velocity profile essen- 
tially shifts downward as LT is reduced. This shift arises, 
in part, as a consequence of differences in the response 
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of the lithosphere to shear stresses acting at the base of 
the region. The net effect of these shear stresses is a 
motion of the lithosphere toward the center of the previ- 
ously glaciated region. (This motion is not evident at all 
locations in Figure 1 because of obscuring contributions 
to tangential motion from other mechanisms.) The thin- 
ner the lithosphere, the less the resistance to the stress 
field, and thus the higher the tangential velocity toward 
the center of the disk which is contributed by the basal 
shear. As vlm is increased to 5 x 10 21 Pa s, the outward 
horizontal speeds in the central region are increased. In 
the far field of the disk center, beyond the main periph- 
eral bulge region, the computed tangential velocity (di- 
rected toward the disk center) increases significantly in 
amplitude as the lower mantle is stiffened. The net effect 
is that relative tangential velocities along the computed 
profiles are particularly sensitive to variations in x >lm- In 
contrast, the relative tangential velocities are insensitive 
to variations in LT. Finally, reducing vvm leads to a sig- 
nificantly broader region in which present-day tangential 
velocities are directed away from the center of the disk. 
Indeed, for the Earth model with uum = 3 x 10 20 Pa s, 
positive tangential velocities extend to a point over 50° 
from the disk center, compared to the corresponding value 
of 16° for the test model. 

Comparing the results of Figures la and lb with those 
obtained using the disk load of maximum radius 8° in- 
dicates that present-day deformation rates of significant 
magnitudes become more concentrated in the central re- 
gion as the length scale of the ice disk is decreased. In 
the case of the test Earth model and the 16° disk load 
the maximum tangential velocities away from and toward 
the disk center are approximately 0.6 and 0.3 mm yr -1 , 
respectively, and the maximum peripheral subsidence is 
1.6 mm yr -1 . The corresponding values for the 8° disk 
load are 1.0, 0.1, and 1.0 mm yr -1 . This result indicates 
that the deglaciation event over Fennoscandia will not 
have a significant impact on the presentrday tangential 
velocity field over, say, North America. 

The computed present-day tangential velocities are less 
sensitive to variations in the lithospheric thickness in the 
case of the 8° disk load than they are in the case of the 
16° disk load. In contrast, the computed radial velocity 
results show greater sensitivity to variations in LT as the 
length scale of the surface mass load is reduced. The 
sensitivity of the computed velocities to variations in vlm 
is much less pronounced for the Fennoscandia scale ice 
sheet since the smaller disk is obviously less efficient in 
driving lower mantle flow than the larger. 

James and Morgan [1990] have computed tangential de- 
formation rates using a fixed-radius disk load model of 
the ancient Laurentide ice sheet, a simplified ocean load 
(extending from 65° to the antipole of the disk), and a 
small number of Earth models. The radial dependencies 
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of the elastic parameters in their Earth models were dis- 
cretized using only five layers, and the model was assumed 
to be elastically incompressible. One consequence of their 
discretization is that the density jump at 670 km depth 
was 1100 kg m -3 , which is almost 3 times the PREM de- 
rived value. The assumption of elastic incompressibility is 
also rather suspect, and, accordingly, we have performed 
calculations, using our Laurentide scale ice disk, which 
explicitly consider this aspect of the James and Morgan 
[1990] models. The solid and dashed lines on Figure 2 
were obtained using the test Earth model; however, in 
the latter case the Earth is assumed to be elastically in- 
compressible. Clearly, the errors incurred by assuming an 
incompressible elastic response, in calculations of present- 
day tangential velocities, can be very large. These errors, 
for the case of the test model, exceed 50% in the central 
region and 35% in the periphery. The errors are smaller in 
the computed peak central uplift rates (about 5%); how- 
ever, they reach 20% in predictions of the peak peripheral 
subsidence. 

James and Morgan [1990] adopted assumptions which 
lead to features in their computed tangential deformation 
rates that are significantly different from those obtained 
in the present analysis. For example, the tangential ve- 
locity profiles computed by James and Morgan [1990, Fig- 
ure 4] do not generally exhibit motion away from the cen- 
ter of the ice disk, within the central region, which is char- 
acteristic of all of the Earth models considered in this sec- 
tion. Indeed, the dominant characteristic of the profiles 
computed by James and Morgan [1990] is a broad, high- 
amplitude, region of negative (i.e., toward the disk center) 
velocity which extends from the load center to a position 
near 70° from this point. The peak amplitude of the neg- 
ative tangential velocities they obtained ranged between 
—2.5 and —4.2 mm yr -1 . (The peak uplift rates for the 
same models were 1. 1-1.5 mm yr -1 .) These tangential 
rates are substantially larger in magnitude than the neg- 
ative velocities presented in this section, regardless of the 
Earth model used, even taking into account the slightly 
higher radial uplift rates. Indeed, our own test calcular 
tions (not shown here) using the ice disk of maximum ra- 
dius 16° and James and Morgan’s [1990] “reference” Earth 
model, predicted a tangential velocity profile with a max- 
imum positive velocity of approximately 1.0 mm yr -1 and 
a peak negative velocity of —0.5 mm yr -1 . We conclude 
that the very large tangential deformation rates directed 
toward the disk center and the absence of central re- 
gions of substantial outward tangential motions, which 
characterize the James and Morgan [1990] results, are ar- 
tifacts of the simplifications adopted in that study. (The 
calculations performed in subsequent sections of this par 
per support this conclusion.) The assumption of a time- 
independent ice disk radius (in addition to that of elastic 
incompressibility) is particularly suspect, and no doubt 
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has a profound influence on the profile of tangential ve- 
locity, especially near the maximum perimeter of the ice 
disk. James and Lambert [1993] have recently suggested 
that the adoption of an incompressible lithosphere in the 
study of James and Morgan [1990] was responsible for the 
prediction of inwardly directed tangential velocities in the 
near field. The results in Figure 2 suggest that other as- 
sumptions adopted by James and Morgan [1990] must also 
contribute significantly in this regard. 

Calculations Using Realistic Ice Models: 

Three Dimensional Deformation Rates 

In Figure 3 we show maps of the predicted present-day 
radial and tangential velocity fields over North Amer- 
ica and Europe generated using the test Earth model. 
The surface mass load incorporated a gravitationally self- 
consistent ocean meltwater component and the ICE-3G 
model for the final late Pleistocene deglaciation event 
[Tushingham and Peltier, 1991], (Independently derived 
maps, based on the same load and Earth model, were 
recently published by James and Lambert [1993].) The 
ICE-3G model has been adapted in order to incorporate 
three full glacial cycles [Mitrovica and Peltier, 1993a]. In 
North America, peak radial velocities are obtained on the 
east coast of James Bay (approximately 13 mm yr _1 ) and 
also near the west coast of Hudson Bay. These locations 
coincide with the two areas in which the (ICE-3G model) 
Laurentide ice sheet thickness exceeded 3 km during the 
last glacial maximum. The main uplift region extends 
over all of Canada and a portion of northeastern United 
States. The secondary uplift peak over the west coast 
of Canada is associated with the disappearance of the 
Cordilleran ice sheet, whereas the deflection of the uplift 
contours toward the north is due to the melting of Arc- 
tic ice. The peak uplift rate over the main Fennoscan- 
dian deglaciation center reaches in excess of 12 mm yr _1 . 
Other regions of uplift are associated with melting over 
Spitzbergen (Svalbard), the Barents and Kara Seas, Ice- 
land, Greenland, and the northern tip of the British Isles. 
Encircling the main regions of uplift are the subsiding 
peripheral bulges. In general, the North American sub- 
sidence rates are somewhat higher than those in Europe, 
as one would expect from the disk load calculations. Sig- 
nificantly, there are several areas in which the subsidence 
rates exceed a few millimeters per year, and these co- 
incide with locations at the periphery of more than one 
ice sheet. Examples are the Davis Strait/Baffin Bay re- 
gion and northwest of Fennoscandia, both of which have 
subsidence rates which reach ~5 mm yr -1 . 

The combination of the ICE-3G ice history and the 
test Earth model has produced peak uplift and subsi- 
dence rates over North America and Fennoscandia which 
are comparable, even though the ice sheet which covered 
the latter was much smaller and melted away much faster 
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than the Lauren tide ice sheet. The reason for this near 
equality is that the characteristic decay time associated 
with the response of the test Earth model increases as the 
spatial scale of the load is reduced (to a limit of about 
700 km). Thus the smaller Fennoscandian ice sheet will 
maintain a relatively greater level of present-day disequi- 
librium. 

Next, consider the tangential velocity fields. Any area 
which is characterized by a local peak in the radial ve- 
locity field, that is a deglaciation center, is subject to 
relatively small tangential deformation rates and encir- 
cled by a region in which the tangential motion is di- 
rected outward (as in Figure 1). In North America the 
pattern is evident near the southern tip of James Bay, 
in Labrador, off the western edge of Hudson Bay, and in 
the Canadian Cordillera. Maximum tangential motions of 
just over 2.0 mm yr -1 are obtained on the northeast coast 
of Labrador, and rates of 1.5 mm yr -1 are evident on Baf- 
fin Island. In Europe this pattern exists over the Gulf of 
Bothnia, the Barents and Kara seas, Iceland, Greenland, 
and the northern reaches of the United Kingdom. The 
peak tangential motions reach ~2 mm yr -1 on the coast 
of Norway and in the Norwegian Sea. Comparable tan- 
gential velocities are also predicted on the western edge of 
the Barents Sea, in the vicinity of Spitzbergen. Encircling 
the whole of Canada, and to the south of the previously 
glaciated areas in Europe, are broad regions in which the 
tangential deformation rates are directed toward the as- 
sociated uplifting regions. These rates are ~0.5 mm yr -1 
in the former region and only slightly less in the latter. 
These amplitudes are significantly larger (approximately 
a factor of two or more) than those recently presented by 
James and Lambert [1993] using an identical Earth model 
along with the ICE-3G deglaciation event. The origin of 
these differences is not known. 

As observed in Figure 1, the melting episode over Lau- 
rentia produces, for predictions based on the test Earth 
model, an appreciable tangential adjustment at very large 
angular distances from Hudson Bay. For the European 
region considered in Figure 3, this contribution should 
be characterized by a roughly uniform tangential mo- 
tion toward the northwest. Thus the tangential mo- 
tions predicted for continental Europe (~0.4 mm yr -1 
to the northwest) arise from the adjustment associated 
with both the local melting events and the melting over 
Laurentia. Indeed, the disk calculations indicate that ve- 
locities of this amplitude could not arise solely from the 
melting of the ice sheets in northern Europe. As a fi- 
nal point in this regard, the Fennoscandian adjustment 
associated with the local melting event will tend to be re- 
inforced, on the northwest side, by the contribution from 
the Laurentide melting episode, and partially cancelled on 
the southeast side for the same reason. As a consequence, 
Norway is experiencing relatively large tangential defor- 
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mation rates of ~2 mm yr -1 , whereas tangential rates 
in Finland and the Baltic Sea (directed to the southeast) 
are about 40% smaller than this value. 

The predicted tangential motions shown in Figure 3 
are a complicated superposition of adjustments associ- 
ated with a number of ice domes. For example, the cor- 
ridor formed by the proglacial system extending from 
Great Bear Lake to Lake Winnipeg is characterized by 
rather weak tangential motions because the eastward and 
westward motions associated with the melting of the 
Cordilleran and Laurentian ice sheets, respectively, al- 
most cancel each other. A similar cancellation is evi- 
dent in Baffin Bay and Davis Strait. In the latter case, 
large gradients (~3 mm yr -1 across a distance of only 
~500 km) occur in the predicted tangential velocity field. 

In the investigations described below we often focus 
on predictions of velocity components along two specific 
great circle profiles; the first extends from Fort Rupert, 
just west of James Bay, to Richmond, Florida; the sec- 
ond extends from Lulea, Sweden to Wettzell, Germany 
(see Figure 3). Fort Rupert and Lulea are near areas 
of peak uplift, while both Richmond and Wettzell are in 
the far field, well south of regions of maximum periph- 
eral subsidence. Figure 4 (solid lines) provides the veloc- 
ity components along these profiles for the calculations 
given in Figure 3. (The values along the Fort Rupert- 
Richmond profile are close to zero because the profile is 
oriented nearly north-south, and the ancient Laurentide 
ice sheet was roughly symmetrical across it. This result 
does not apply to the European profile since that profile 
is neither an approximate symmetry axis for the uplifting 
Fennoscandian region, nor is it oriented north-south.) 

We now investigate the sensitivity of the results in Fig- 
ures 3 and 4 (solid lines) to variations in both the surface 
mass load and the Earth model. Figure 5 shows predic- 
tions (analogous to those presented on Figure 3) gener- 
ated using an ice load adapted, as described above, from 
the ICE-1 deglaciation history [Peltier and Andrews, 1976). 
Comparing Figures 3 and 5, we see a marked similarity 
in the general form of the radial and tangential veloc- 
ity fields predicted on the basis of the two ice models. 
(As we have observed the same is not necessarily true of 
predictions based on Earth models with different mantle 
viscosity profiles.) Important differences are, however, 
apparent. For example, the peak uplift rates over Hud- 
son Bay and Fennoscandia are larger (by about 20%) for 
predictions based on the ICE-1 model. In the case of the 
ICE-1 results the region of peak uplift in Canada extends 
east-west over a broad area which encompasses southern 
Hudson Bay (and which is somewhat north of the peak 
uplift, at Fort Rupert, on Figure 3). Furthermore, in 
Fennoscandia, the major axis of the uplift field predicted 
using the ICE-1 model is, in contrast to Figure 3, aligned 
roughly parallel to the Baltic Sea and the Lulea- Wettzell 
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profile. These differences, due to variations in the ice 
sheet geometry, are also evident on the tangential veloc- 
ity maps. In the case of predictions generated using the 
ICE-1 model the central region in which tangential veloc- 
ities are very small (and from which the velocity vectors 
emanate outward) is located slightly northward (~300 
km) in the Hudson Bay region and slightly southward 
(~200 km) in the Baltic Sea area of the corresponding lo- 
cations for the ICE-3G predictions. Differences between 
Figures 3 and 5 are also significant over Greenland, where 
the ICE-1 model adopts a relatively insignificant deglaciar 
tion event, and over the Arctic region north of Europe. 
In the latter area the ICE-3G model includes a Barents 
Sea ice complex, whereas the ICE-1 history does not. 

The velocity components along the Fort Rupert to 
Richmond and Lulea to Wettzell profiles, generated using 
the ICE-1 and ICE-3G models (Figure 4), reflect the vari- 
ations described above. (Notice, for example, the relative 
shift between the Ug component of velocity computed us- 
■> ing the two ice models; the slight oscillation evident in 
the velocity profiles computed using the ICE-1 model is 
due to the very coarse spatial discretization which char- 
acterizes that model; see Peltier and Andrews [1976].) The 
results indicate that predictions of three-dimensional ve- 
locities at specific geographic sites are sensitive to varia- 
tions in the ice mass history. This sensitivity may not be 
surprising since we have already observed (in the context 
of Figure 3) that velocity predictions at a particular site 
are strongly dependent on the location of the site relative 
to the spatial geometry of the ice load. (The location of 
the two great circle profiles was based on an inspection 
of the ICE-3G based results in Figure 3.) Clearly, analy- 
ses based on a comparison of predicted three-dimensional 
velocities with observational constraints (for example, to 
infer mantle rheology) must recognize this sensitivity. 

In Figure 6 we investigate the sensitivity of the predic- 
tions in Figure 3 to errors in the ocean loading model. 
The solid lines denote velocity components generated 
using a gravitationally self-consistent ocean load. The 
dashed and dotted lines represent results generated us- 
ing simplified models of the ocean mass redistribution. 
The dashed line is computed by adopting the so-called 
eustatic approximation. That is, it is assumed that the 
ocean mass load, which complements the mass fluctuation 
in the late Pleistocene ice sheets, is independent of po- 
sition [e.g., Mitrovica and Peltier, 1991]. The dotted lines 
represent results generated by assuming no ocean load 
component. 

We conclude that radial velocities can be computed ac- 
curately using the eustatic approximation, except in re- 
gions within, or in the vicinity of, oceans. The largest er- 
rors in the Fort Rupert-Richmond profile (~0.5 mm yr _1 ) 
occur at the edges of the profile, that is, in the vicinity 
of Hudson Bay and the Atlantic Ocean. Even more pro- 
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nounced errors (~2.0 mm yr -1 ) are evident on the first 
1000 km of the profile between Lulea and Wettzell, which 
skirts and sometimes enters the western coasts of the Gulf 
of Bothnia and the Baltic Sea. These errors are due to 
the neglect of hydroisostatic effects. In a gravitationally 
self-consistent calculation of the ocean redistribution the 
on-going removal of ocean mass in the region of post- 
glacial uplift, as a consequence of this uplift, represents 
a continuously increasing negative mass load on the sys- 
tem. This load will tend to increase the present-day ra- 
dial velocity in the region and its vicinity. The eustatic 
approximation does not incorporate this effect since, in 
this case, the ocean load at any particular time represents 
some globally averaged distribution of the complement of 
the total ice mass. Indeed, since the ICE-3G deglaciar 
tion event ends at 5 kyr B.P., the eustatic approximation 
implies no change in the ocean load subsequent to this 
time. 

The influence of hydroisostasy on the computed pre- 
sent-day tangential velocities is not insignificant. Along 
the profiles the amplitude of the error which results from 
using the eustatic approximation can be as large as 0.1- 
0.2 mm yr -1 , with the peak obtained at the point of maxi- 
mum tangential speed on the profiles (~1 mm yr -1 ). The 
sign of the error is opposite for the two profiles. For the 
Fort Rupert to Richmond profiles, hydroisostatic effects 
diminish the tangential speed by as much as 0.1 mm yr -1 
in this region (~500 km south of Fort Rupert), while the 
same effects add approximately 0.2 mm yr -1 to the tan- 
gential speed of the 750 km point on the Lulea- Wettzell 
profile. The reason for this difference is that the posi- 
tion of the ocean load relative to the profile is different 
in the two cases. The 500-km point on the Fort Rupert 
to Richmond profile is well outside the Hudson Bay area. 
Therefore a negative mass load on Hudson Bay will tend 
to produce (using the test Earth model) an incremental 
tangential motion toward Hudson Bay. The 750-km point 
on the Lulea to Wettzell profile is, in contrast, essentially 
in the very midst of the ocean load. Thus the negative 
mass load will produce incremental tangential velocities 
away from the region. 

Assumptions made concerning the modeling of the 
ocean load are important when considering tangential 
velocities associated with specific coastal sites (see Fig- 
ure 7). At some coasted sites the eustatic approximation 
is accurate (e.g., sites 5-7, on the east coast of the United 
States), whereas at some others (e.g., site 16, Helsinki) 
it is no more accurate than the assumption of no ocean 
load. The errors incurred using either a nonexistent or 
gravitationally inconsistent ocean load component are ev- 
ident in both the amplitude and direction of the predicted 
tangential motions. 

In Figure 8 we consider results generated by varying 
the lithospheric thickness of the test model. The con- 
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elusions based upon the simple disk load calculations are 
valid for the interpretation of the sensitivities evident on 
Figure 8. For example, the computed radial velocity pro- 
files in Fennoscandia are more sensitive to variations in 
LT than are the corresponding North American profiles. 
The effect of increasing LT on the computed radial veloc- 
ity profiles is to reduce the present-day peak uplift and 
subsidence rates. Indeed, the peak relative radial veloc- 
ities on the European and North American profiles are 
reduced, respectively, by 3.5 and 1.5 mm yr -1 , as LT is 
increased from 96 to 145 km. Increasing the lithospheric 
thickness tends to smooth the three-dimensional velocity 
profiles. Increasing LT also results in a relatively uniform 
positive perturbation in the tangential velocity field. Rel- 
ative tangential velocities are clearly not strongly affected 
by variations in LT. 

In Figures 9 and 10 we consider the sensitivity of 
the three-dimensional velocity predictions to variations 
of vlm and uum, respectively. As in Figure 1, increasing 
ulm leads to an increase in the peak present-day uplift 
rate and an outward migration of the position of max- 
imum subsidence on the peripheral bulge. Decreasing 
uum results in predictions with lower peak uplift and 
subsidence rates. Since the ice sheet which existed over 
Laurentia was much larger than that which existed over 
Fennoscandia, the computed radial uplift rates along the 
Fort Rupert-Richmond profile tend to be more sensitive 
to variations in u^m and less sensitive to variations in 
than are the uplift rates for Lulea-Wettzell. 

It is difficult to compare directly the tangential veloc- 
ity profiles of Figures 9 and 10 with corresponding results 
obtained using simplified disk load models. The realistic 
ICE-3G ice sheets exhibit profound asymmetries, both 
on local and global length scales. As a result, the specific 
profiles we are considering may be affected in a compli- 
cated manner; certainly, the profiles are not constrained 
to have zero tangential velocities at their starting points. 
Furthermore, the influence of asymmetries in the surface 
mass load on the computed profiles will be model depen- 
dent since tangential motions are sensitive to variations in 
the Earth model characteristics, even at great distances 
from the load. Similarities between the tangential veloc- 
ity predictions shown in Figures 9 and 10, and the cor- 
responding predictions from disk load models, are more 
obvious when one considers relative tangential velocities. 

Figure 11 shows the present-day values for Ug along 
each of the Fort Rupert-Richmond and Lulea-Wettzell 
profiles relative to the rate at the starting point of each, 
for the Earth models considered in Figure 9. Increas- 
ing v lm produces an increase in the (relative) outward 
tangential velocity in the near field and, in the case of 
the larger ice sheet, a significant increase in the (relative) 
tangential velocity toward the previously glaciated area 
in the far field. (In the case of ulm = 10 21 Pa s, Rich- 
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mond is, in contrast to the test Earth model results, pre- 
dicted to be moving away from the previously glaciated 
region.) The relative tangential velocities in the far field 
for the European profiles exhibit only weak sensitivity to 
variations in vlm, whereas the absolute velocities (Fig- 
ure 9b) are much more sensitive to these variations. This 
difference is due to the adjustment in Fennoscandia, asso- 
ciated with the melting of the Laurentide ice sheet, being 
sensitive to variations in 

The adjustment due to Laurentide deglaciation must 
also be influencing the tangential velocity results in Fig- 
ure 10b, in which the sensitivity to variations in vum of 
the computed profiles is explored. However, this influence 
is not as apparent in Figure 10b, as it was in Figure 9b, 
mainly because the adjustment associated with the melt- 
ing of the (local) Fennoscandian ice complex is also sen- 
sitive to variations in i/y M . Relative tangential velocities 
along the profiles of Figure 10b arise mainly from local 
melting events. Once again, an important observation is 
that a reduction in vum leads to a significant enlarge- 
ment of the region characterized by tangential motions 
directed outward toward the far field. 

To gain insight into the sensitivity of a particular ob- 
servable to the detailed radial variation of rheology re- 
quires the calculation of Frechet kernels (equation (2)). 
Frechet kernels for predictions of present-day three-dimen- 
sional deformation rates at nine different sites on both the 
Fort Rupert-Richmond and Lulea-Wettzell profiles are 
shown on Figure 12. The kernels were computed using 
the test Earth model, the ICE-3G adapted ice history, 
and a gravitationally self-consistent ocean load. 

Consider the kernels for site A, situated at Fort Rupert, 
which is near the location of maximum uplift rate. The 
kernel for the radial velocity indicates that the velocity 
has a nonnegligible sensitivity to variations in viscosity at 
depths extending down to the CMB, although the max- 
imum sensitivity is clearly located in the top half of the 
lower mantle and in the upper mantle. This suggests that 
the sensitivity of the predictions of U r at Fort Rupert to 
variations in vlm (Figure 9a) should be ascribed mainly 
to a sensitivity to variations in rheology above 1600 km 
depth. In contrast, the sensitivity of the same radial ve- 
locity prediction to variations in vum (see Figure 10a) 
reflects a sensitivity which is more or less constant over 
the entire upper mantle region. The ratio of the area un- 
der the radial velocity kernel for Fort Rupert in the lower 
mantle to the corresponding area in the upper mantle 
is approximately 2:1. This ratio remains relatively con- 
stant for the sites B and C, which are also situated in 
the central uplift region of Laurentia. Indeed, the main 
difference in the radial velocity kernels as one moves from 
site A to site C is a progressive and rather uniform re- 
duction in amplitude. This indicates that radial velocity 
data obtained at such sites will not provide independent 


14 



information regarding the mantle viscosity profile. 

The U e kernels indicate that the predicted Ug values 
for sites A-C are much less sensitive to variations in the 
rheology of the lower mantle than to variations in the rhe- 
ology of the upper mantle. (See also Figures 9 and 10.) 
Indeed, the sensitivity to variations in viscosity within 
the lower mantle is limited to the very shallowest por- 
tions of this region. Furthermore, a comparison of these 
kernels with those for the radial velocity predictions indi- 
cates that uplift rate data from central regions will pro- 
vide substantially more information regarding deep man- 
tle rheology than tangential velocity data from the same 
area. 

The Frechet kernels in Figure 12, sites D-I, suggest 
that the detailed depth-dependent sensitivity of the three- 
dimensional velocity components varies significantly as 
one moves from the edge of the uplift region, through 
the peripheral bulge, and on to the far field. Three- 
dimensional velocity data obtained from each of these 
regions will provide independent constraints on the ra r 
dial profile of mantle viscosity. 

As one proceeds from a site just inside the edge of the 
uplift region (site C), to one just outside the edge (site D), 
the radial velocity kernels continue to trend toward less 
positive values. This trend is most evident in the shal- 
lowest regions of the mantle, and it leads, for site D, to 
negative values of the kernel within the upper mantle. 
This suggests that at some point between the sites C and 
D, near the very edge of the predicted uplift region, the 
area under the associated U r kernel in the upper mantle 
will vanish. Thus the radial velocity computed at this 
point using the test Earth model should be insensitive 
to variations in v\jm- This location coincides with the 
“hinge point,” 900 km from Fort Rupert, in the bottom 
frame of Figure 10a. The radial velocity kernel in Fig- 
ure 12, site D, indicates that the U r prediction at site 
D using the test Earth model will have a nearly uniform 
sensitivity to variations in viscosity anywhere in the lower 
mantle. Thus the migration of the perimeter of the uplift 
region, observed in numerical studies of glacial isostatic 
adjustment in North America in which vlm is varied [e.g., 
Cathles, 1975; Peltier and Andrews, 1976], is sensitive to 
variations in the rheology of the entire lower mantle. 

Continuing from outside the edge of the uplift region 
(site D) to a point somewhat beyond the location of max- 
imum peripheral subsidence rate (site E), the emergence 
of negative values on the radial velocity kernels contin- 
ues and extends through the top half of the lower mantle. 
Within this range of sites the prediction of radial velocity 
becomes increasingly sensitive to variations in the upper 
mantle viscosity. Between sites D and E the Frechet ker- 
nel in the lower mantle evolves from being nearly uniform 
and positive to being mostly negative and peaking at the 
top of the region. At some point between sites D and E, 
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but closer to site E and, very close to the position of max- 
imum subsidence on the peripheral bulge, the area under 
the radial velocity kernels in the lower mantle vanishes. 
As one proceeds from site E through to the outer flank of 
the peripheral bulge (site F) and then on to the far-field 
(sites G-I), the region exhibiting the peak sensitivity to 
variations in rheology migrates to greater depths in the 
mantle. In the outer regions of the peripheral bulge (sites 
E to F) the predicted radial velocity is similarly sensitive 
to variations in the viscosity of the top half of the lower 
mantle and in the viscosity of the upper mantle. In the 
far field, however, the predictions are almost totally in- 
sensitive to variations in viscosity anywhere in the upper 
mantle. The sensitivity of the radial velocity prediction 
at Richmond to variations in I'm, evident in the bot- 
tom frame of Figure 9a, may be ascribed to a sensitivity 
which peaks in the lower mantle between depths of 800 
and 2200 km. 

The change in the radial velocity kernels between sites 
D and I provides a detailed picture of the sensitivity of 
the dynamics of the peripheral bulge to variations in lower 
mantle rheology. The predicted migration of the periph- 
eral bulge with variations in vlm is a consequence of 
the reversal in the area under the radial velocity kernels, 
within the lower mantle, as one proceeds from sites on the 
near field side of the point of maximum peripheral subsi- 
dence (site D) to sites on the far-field side of that point 
(sites E-G). The reversal is associated with significant 
changes in the sensitivity of the predicted radial velocity 
to variations of rheology at the top of the lower mantle, 
as one moves along the profile. This evolving sensitivity 
produces a “twisting” of the predicted profile, within the 
periphery, as is varied (clockwise for positive variar 
tions) and, consequently, the migration of the peripheral 
bulge. This migration is evident in the results described 
herein (Figure 9a) and in a number of classic analyses 
[e.g., Cathles, 1975; Peltier and Andrews, 1976J. 

The evolution of the FVechet kernels for Ug along the 
North American profile is less complicated than the evolu- 
tion in the radial velocity kernels. As one proceeds from 
sites B and C, which straddle the point on the profile 
characterized by the maximum positive tangential veloc- 
ity, through the peripheral bulge and the far field, the 
peak amplitude of the kernel diminishes at the same time 
as the deep mantle sensitivity of the tangential velocity 
increases. For example, predictions at sites within the 
main peripheral bulge region (D-F) are sensitive to vari- 
ations in the lower mantle rheology only in the top 500- 
700 km or so of that region, whereas the associated calcu- 
lations at far-field sites (G-I) are sensitive to variations 
in viscosity throughout the lower mantle. This trend is 
evident on the middle frame of Figure 9a, where the dis- 
crepancy between the three profiles generally increases as 
one moves toward Richmond in the far field. One of the 
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main results of the forward calculations described above 
was the prediction of extensive regions in the far field of 
the Laurent ide ice sheet which are characterized by mo- 
tions toward the previously glaciated region. The striking 
sensitivity of this computed motion to variations in vlm 
can therefore be ascribed to a sensitivity which persists 
over the entire lower mantle. Indeed, an examination of 
the kernel for the Ug prediction for Richmond indicates 
that variations in the rheology of the top half of the man- 
tle are only a factor of 2 more efficient at perturbing the 
tangential velocity computed using the test Earth model 
than variations in the viscosity in the bottom half of the 
mantle. 

The Frechet kernels for the tangential velocity compo- 
nents in Figure 12, sites A-I, are characterized by uni- 
formly negative values in the upper mantle. This negar 
tivity indicates that a reduction in vum should produce a 
positive shift in the tangential velocity components pre- 
dicted for the Fort Rupert-Richmond profile (or a broad- 
ening of the central zone of outward tangential motions, 
Figure 10a) relative to predictions generated using the 
test Earth model. This shift will be particularly signif- 
icant for the case of the Ug component, and it will be 
most pronounced in the near field, where the amplitudes 
of the kernels within the upper mantle peak (sites A-F). 
Figure 12 (sites A-F) indicates that the broadening is 
particularly sensitive to variations in rheology in progres- 
sively shallower regions of the upper mantle. 

The evolution of the Frechet kernels in Figure 12, sites 
J-R, is reminiscent of the trends evident along the North 
American profile; however, the kernels for the European 
sites tend to be “scaled” toward the surface in compari- 
son to the corresponding kernels for North America. This 
result is not surprising given the difference in the spa- 
tial scales of the Laurentian and Fennoscandian ice com- 
plexes, and it has several interesting implications. First, 
in contrast to the North American edge site D, the radial 
velocity kernel for the edge site M indicates a sensitiv- 
ity to lower mantle viscosity which is mainly ascribable 
to the top half of that region. Second, a migration to 
greater depths of the peak value of the radial velocity 
kernels as one proceeds through the far field, evident in 
the North American results, is not nearly as pronounced 
in the European results. The amplitude of the radial ve- 
locity kernel for Wettzell (Figure 12, site R) is largest over 
depths in the range 400-1300 km, in contrast to the cor- 
responding range of 800-2200 km for the Richmond site. 
Third, as one proceeds from the central to far-field sites, 
the region of peak sensitivity exhibited by the kernels for 
the tangential velocity Ug in Europe never extends below 
the top 300 km of the lower mantle. 

The tangential motions predicted for the North Ameri- 
can profile tend to be more sensitive to variations in vum 
than the corresponding predictions over Europe (see Fig- 
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ure 10). The kernels shown in Figure 12 indicate that 
this difference is due not only to the difference in the 
peak sensitivities evident on the kernels for sites at cor- 
responding points on the two profiles (sites D and M, 
for example), but also to the different depth-dependent 
sensitivities which characterize these kernels within the 
upper mantle. 

The FVechet kernels for tangential velocity in Europe 
exhibit nonzero sensitivity to variations in the viscosity 
in the deep mantle extending to the CMB. This is a con- 
sequence of the contribution to the predicted tangential 
velocities in Europe from the melting of the Laurentide 
ice sheet over Canada. The deep mantle sensitivity ev- 
ident in Fennoscandia should be partitioned into 0 and 
ip components since the great circle arc extending from 
Fennoscandia to Laurentia is roughly in a northwest di- 
rection, viewed from Fennoscandia. The results of Fig- 
ure 12, sites J-R, reflect this partitioning. 

Calculations Using Realistic Ice Models: 

Comparison with VLB I Results 

Ryan el aL [1993] have analyzed Mark III geodetic VLBI 
data obtained from 1979 to 1991 to determine the evolu- 
tions of a large number of baselines. Their “global” so- 
lution (GLB868) included 864,359 group delay measure- 
ments obtained in 1648 observing sessions. In this section 
we consider predictions of velocity vectors (equation (1)) 
generated using a suite of Earth models (see Table 1) for 
a subset of the baselines included in the GLB868 solu- 
tion. The subset is chosen to satisfy the following crite- 
ria: (1) At least one of the sites defining each baseline is 
situated within the region enclosed by the outer edge of 
the peripheral bulge associated with the melting of either 
the Laurentian or Fennoscandian ice complexes. (2) The 
sites defining the baseline are both located on the North 
American plate or, alternatively, both on the European 
plate. This restriction ensures that we do not have to cor- 
rect the baseline velocity results for relative plate motions 
or to consider the errors in this correction. (3) The base- 
lines do not include sites (e.g., Penticton) near the North 
American-Pacific plate boundary, which may be subject 
to significant deformations. These criteria are satisfied by 
28 baselines, 23 of which are located in North America. 
The location of the baselines are given in Figure 13. 

Tushingham [1991] has computed baseline rates associ- 
ated with glacial isostatic adjustment under the assump- 
tion that the sites defining the baseline experience only 
radial motions. In Table 2 we present calculations us- 
ing the ICE-3G adapted ice history and a gravitationally 
self-consistent ocean load intended to assess the validity 
of this assumption (see also James and Lambert [1993] and 
Mitrovica et aL [1993]). The format of Table 2a is as 
follows: The columns are separated into three triplets, 
each of which refers to results obtained with a different 
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Earth model. The three columns associated with each 
Earth model, labelled TOT, RAD and TAN, represent 
the total predicted baseline velocity component and its 
contribution from radial and tangential site velocities, re- 
spectively. Table 2a includes two representative baselines, 
Algonquin Park-Gilmore Creek and OnsalarWettzell. In 
both cases, three rows corresponding to the L , T, and 
V components of the predicted baseline velocity compo- 
nents are included (as labeled). Table 2b provides the 
three components of present-day site velocity for the four 
sites represented on the two baselines, for each of the 
Earth models. 

The analysis of Tushingham [1991] was limited to base- 
lines within North America and (radial) velocity predic- 
tions based on the test Earth model only. In the case 
of the test model results the tangential motions con- 
tribute 15% of the predicted length component of the Al- 
gonquin Park-Gilmore Creek baseline velocity. However, 
the results in Table 2a indicate that the Algonquin Park- 
Gilmore Creek baseline length rates are not dominated by 
radial site motions for Earth models other than the test 
model. Furthermore, the rate of change of the Onsala- 
Wettzell baseline length is dominated by the contribution 
from tangential motions regardless of the Earth model 
considered. (The results of James and Lambert [1993] 
also indicate the importance of rebound-related tangen- 
tial motions to baseline vector velocity.) The Tushingham 
[1991] approximation is therefore not generally valid. 

The relatively small contribution from tangential mo- 
tions to the length component of the Algonquin Park- 
Gilmore Creek baseline velocity, in the case of the test 
Earth model, is a result of a relatively rare match be- 
tween the tangential velocities along the direction of the 
baseline at the two end sites. In this case, both sites 
are moving in the same direction along the baseline, with 
approximately the same velocity (see Figure 3). For com- 
putations using model 1, which has a weak upper mantle 
viscosity, these velocities were both significant and di- 
rected oppositely (see Table 2b). For model 5 which has 
a stiff lower mantle, the velocity components were in the 
same direction but of very different magnitude. 

Tangential velocities always contribute a dominant por- 
tion of the predicted length rate of the OnsalarWettzell 
baseline for two reasons. First, since the baseline is rel- 
atively short (920 km), radial motions at the end sites 
project weakly onto the baseline direction. Second, both 
Onsala and Wettzell are located in regions which are not 
characterized by large radial motions (Figure 13). This 
is in contrast, for example, to the Algonquin Park site. 

The results summarized in Table 2a have important 
implications for the sensitivity of the baseline velocity 
components to variations in rheology. For regional base- 
lines the vertical component of baseline velocity tends to 
be dominated by the radial motions of the end sites. Thus 
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the sensitivity of this component to variations in mantle 
viscosity at depth should be approximately proportional 
to the difference between the Frechet kernels for the rar 
dial velocity at the two sites. The FYechet kernel for the 
transverse component of baseline velocity is a weighted 
sum of the Frechet kernels for the tangential components 
of velocity at each of the two sites. Finally, the Frechet 
kernels for the length component of the baseline velocity 
will incorporate (i.e., “mix”) the sensitivities associated 
with both radial and tangential motions at the two sites. 
The results in Table 2a suggest that this mix will gen- 
erally be dominated by the sensitivity associated with 
tangential motions for the European baselines. The same 
is true, although to a lesser extent, for the Algonquin 
Park-Gilmore Creek baseline. 

To illustrate these points, in Figure 14 we show FYechet 
kernels for a characteristic subset of the North American 
baselines and for all of the European baselines. The ker- 
nels for the vertical component of baseline velocity indi- 
cate that this component has a depth-dependent rheolog- 
ical sensitivity reflective of the sensitivities inherent to 
the radial motions associated with the two sites defining 
the baseline. For example, Algonquin Park and Gilmore 
Creek are positioned with respect to the velocity field 
at sites corresponding closely to sites C and G on Fig- 
ure 9. The difference between the radial velocity kernels 
on Figure 12, sites C and G, yields a kernel with essen- 
tially the same form as the dashed line on Figure 14 for 
the Algonquin Park-Gilmore Creek baseline. The Frechet 
kernel for the vertical component of the Algonquin Park- 
Gilmore Creek baseline velocity indicates that this com- 
ponent will be about twice as sensitive to a variation in 
V LM from the test model value than to an analogous vari- 
ation in v\jm (see also Table 2a). Furthermore, the ker- 
nel indicates a sensitivity which becomes progressively 
stronger as one proceeds from deeper to shallower regions 
of the lower mantle. 

The Frechet kernel for the length component of the Al- 
gonquin Park-Gilmore Creek baseline velocity indicates 
that this component is progressively more sensitive to 
variations in rheology as one proceeds from deeper to 
shallower regions of the lower mantle and relatively in- 
sensitive to variations in the rheology of the lower man- 
tle. A comparison of this kernel with those in Figure 12, 
sites C and G, indicates that the sensitivity of the base- 
line length rate is associated mainly with the sensitivity 
inherent to the tangential motions at the end sites. The 
same conclusion generally applies to the other baselines 
considered in Figure 14. 

It is difficult to identify universally valid trends in the 
FYechet kernels shown in Figure 14. The baseline com- 
ponents, though not necessarily insensitive to variations 
in the rheology of the test Earth model in the deepest 
regions of the mantle, are clearly more sensitive to such 
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variations in the top half of the mantle. Perhaps the most 
significant aspect of the results in Figure 14 is that the set 
of baseline data exhibit significant differences in their sen- 
sitivity to depth-dependent variations in rheology. These 
differences indicate that the ensemble of baseline veloc- 
ity data, given sufficiently accurate information regarding 
the ice loading history, has the potential to provide signif- 
icant constraints on the rheology of various regions within 
the upper mantle and the top half of the lower mantle. 

Table 3 provides the L, T, and V components of 
the baseline evolution computed using the Earth models 
listed in Table 1, the ICE-3G adapted ice history and a 
gravitationally self-consistent ocean load. (Independently 
derived predictions, for 10 baselines, for the specific case 
of the test Earth model, have also recently been published 
by James and Lambert [1993].) It is clear from both Table 3 
and Figure 14 that the predicted baseline velocity compo- 
nents are a strong function of the radial viscosity profile 
of the Earth model adopted in the calculations. The re- 
sults of the last section indicate that these velocities will 
also be sensitive to variations in the adopted ice model. 
Accordingly, it may be possible to discriminate between 
specific combinations of ice history and Earth rheology on 
the basis of the fit between the predicted baseline veloc- 
ity components and the components estimated from the 
analysis of the VLBI data set. In the remainder of this 
section we perform such an analysis using predictions of 
the length rate of the five European baselines (Table 3). 
An analogous study which focuses on the North American 
baselines is given by Mitrovica el aL [1993]. 

Figure 15 is a plot of the x 2 residual per degree-of- 
freedom (henceforth “reduced x 2 ”) of the VLBI deter- 
mined length rates for the European baselines relative to 
the predictions obtained using the indicated Earth model. 
(The uncertainties of the baseline rate estimates used for 
the x 2 calculations and presented below are those given 
by Ryan el aL, [1993], i.e., the statistical standard devia- 
tions of the rates scaled by the square-root of the reduced 
X 2 residual of the series of baseline length estimates from 
the best-fit line.) The reduced x 2 computed assuming 
that there is no influence on baseline rates due to glacial 
isostatic adjustment (henceforth “reduced x 2 of the null 
hypothesis”) is ~2.4. The dominant contributor to the 
reduced x 2 for the null hypothesis is the length rate for 
the rather precisely determined (—0.6 ± 0.2 mm yr -1 ) 
Onsala-Wettzell baseline. The Onsala-Medicina baseline, 
with an observed shortening of 2.1 ± 1.3 mm yr -1 , also 
contributes importantly to the x 2 residual for the null hy- 
pothesis. Clearly, corrections to the VLBI observations of 
baseline length changes in Europe, for the predicted sig- 
nature of glacial isostatic adjustment, will lower the x 2 
residual only if these predictions account for some of the 
shortening evident on the Onsala-Wettzell baseline and, 
to a lesser extent, on the Onsala-Medicina baseline. It is 
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clear from Figure 15 (and also Table 3) that predictions 
based on a number of Earth models, together with the 
ICE-3G adapted ice loading history, are successful in this 
respect. 

Model 2, for example, which is characterized by LT = 
120 k = 5 x 10 20 Pa s, and i 'lm = 2x 10 21 Pa s, 
yields (in combination with the ICE-3G model) predic- 
tions of length changes for both the Onsala-Wettzell 
and Onsala-Medicina baselines of —0.7 mm yr -1 . If we 
use these values to “correct” the observed VLBI base- 
line length rates, we obtain residual values of —0.1 ± 
0.2 mm yr -1 for Onsala-Wettzell and — 1.4±1.3 mm yr -1 
for Onsala-Medicina. Neither residual is significantly dif- 
ferent from zero; the associated reduced x 2 f° r the Eu- 
ropean baseline length changes is only 0.4. It is perhaps 
significant that Model 2 is consistent with the inference 
of Mitrovica and Peltier [1993b] based on an inversion of 
the Fennoscandian relaxation spectrum. These authors 
found, in the case of LT = 120 km, that their preferred 
two-layer model consisted of i/ UM = 4 x 10 20 Pa s and 
vlm = 2x 10 21 Pa s. 

Model 2 does not, however, uniquely “fit” the VLBI 
estimates. Within the limited class of models considered 
in Figure 15 there are several which yield predictions of 
baseline length changes in Europe which are consistent 
with the VLBI values. Indeed, reduced y 2 values below 
1.0 are obtained for models derived from the test model 
either by reducing ulm to 10 21 Pa s or below, by reducing 
uum to a value between 3 and 7 x 10 20 Pas, or by increas- 
ing LT above about 150 km. Of course, the nonuniqueness 
is even more severe when one takes into account the un- 
certainty in the ice history. A formal geophysical inver- 
sion is required to identify the non-uniqueness inherent 
to the inference; however the results in Figure 15 rep- 
resent the solution of the more basic existence problem. 
That is, we have identified at least one Earth model/ice 
model combination and, in fact, several such combina- 
tions, which fit to the VLBLdetermined baseline length 
rates for the five European baselines. Given that none 
of the VLBI baseline rate estimates represent a “detec- 
tion” at the level of greater than 3 -o (and most are less 
than 2 -a), a more rigorous inversion, or a more labored 
interpretation, would not be warranted. 


Summary and Final Remarks 

We have applied the spectral formalism outlined by 
Mitrovica et aL, [this issue] to predict present-day three- 
dimensioned deformation rates and baseline rates driven 
by the ice and ocean loads associated with the late Pleis- 
tocene glacial cycles. One of the mean conclusions of the 
anedysis is that the predictions are sensitive to variations 
in both the viscosity profile of the Earth model and the 
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details of the late Pleistocene deglaciation event. This 
sensitivity suggests that “standard” corrections for the 
signal associated with glacial isostatic adjustment should 
not be applied to geodetic data in an attempt to remove 
the signal. Rather, one must consider the potentially 
large range of signals associated with the phenomenon 
(see, e.g., Table 3), or simultaneously estimate the glacial 
isostatic adjustment signal. If the latter is possible then 
the sensitivities indicated by the present results suggest 
that attempts may be made to refine constraints on ei- 
ther the ice load history or the Earth rheology, or, more 
likely, combinations of the two. 

James and Lambert [1993] have recently stated that pre- 
dictions of tangential motions are much more sensitive 
to details of the viscoelastic Earth models than are rar 
dial motions. This generalization is not supported by the 
results of the present study, which suggest that the sen- 
sitivities are a strong function of the geographic location 
of the site and the specific parameter of the Earth model 
being varied. Consider, for example, Figure 9, site A. In- 
creasing uim by a factor of 5 alters the predicted radial 
velocity by 9.7 mm yr -1 and the horizontal velocity by 
only 0.6 mm yr^ 1 . 

We did not consider transcontinental baselines in the 
present analysis since those baselines will also be sensitive 
to plate motions. The contributions to the velocities of 
sites in North America and Europe due to glacial isostatic 
adjustment are nonnegligible, and should be taken into 
account in an analysis of the time behavior of transcon- 
tinental baselines. It is important to point out, however, 
that the continuity of the tangential deformation field as- 
sociated with melting of Laurentide scale ice sheets was 
based on spherically symmetric Earth models. Disconti- 
nuities in lithospheric strength associated with, e.g., plate 
boundaries, may well influence this trend. 
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Figure 1. Present-day deformation rates computed using 
the ice disk models with maximum radius (a) and (b) 16° 
and (c) and (d) 8°. The solid line on each frame is gen- 
erated using the test Earth model. The remaining lines 
are generated using an Earth model constructed by varying 
a single parameter of the test model: LT=96 km (dashed 
line); wjm — 3 x 10 20 Pa s (dashed-dotted line); and 
vlm — 5 x 10 21 Pa s (dotted line). As a consequence of 
symmetry, the tangential velocity is identically zero at the 
center of the disk location and at the antipole. A positive 
tangential velocity refers to motion directed away from the 
center of the disk load. 


Figure 2. As in Figures la and lb, except that the dashed 
line was computed using the test Earth model with an elasti- 
cally incompressible rheology. A positive tangential velocity 
refers to motion directed away from the center of the disk 
load. 


Figure 3. Plots of the predicted (bottom) present-day ra- 
dial and (top) tangential velocity over North America and 
Europe computed using the test Earth model and a sur- 
face mass load composed of the ICE-3G adapted ice history 
and a gravitationally self-consistent ocean load. The solid 
dots (top frames) refer to the sites: Gilmore Creek, Alaska 
(G); Fort Rupert, Quebec (F); Algonquin Park, Ontario (A); 
Richmond, Florida (R); Onsala, Sweden (0); Lulea, Sweden 
(L); and Wettzell, Germany (W). 


Figure 4. (a) The present-day three-dimensional veloci- 

ties computed along the profile from Fort Rupert to Rich- 
mond using the test Earth model and a gravitationally self- 
consistent ocean mass redistribution. The solid and dashed 
lines were computed using ice loads adapted from, respec- 
tively, the ICB-3G and ICE-1 deglaciation histories, (b) As 
in Figure 4a, except for the profile extending from Lulea to 
Wettzell. 


Figure 5. As in Figure 3, except for an ice load adapted from 
the ICE-1 deglaciation history. 
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Figure 6, As in Figure 4, except all calculations use an 
ice load adapted from the ICE-3G deglaciation history and 
the test Earth model. The solid, dashed and dotted lines 
were generated by assuming a gravitationally self-consistent 
ocean mass redistribution, a eustatic ocean mass redistribu- 
tion, and no ocean load component, respectively. 


Figure 7. Vector plots of present-day tangential velocities at 
selected coastal sites in (a) eastern North America and (b) 
Europe, computed using the test Earth model and the ICE- 
3G adapted ice loading history. The vectors whose heads are 
given by solid triangles, open triangles, or two lines, refer to 
results generated using a gravitationally self-consistent, eu- 
static, or no ocean load. The numbers refer to the sites: 
1, Fort George, Quebec; 2, Indian Harbour, Newfound- 
land; 3, Die d’Anticosti, Quebec; 4, Quebec City, Quebec; 
5, Boston, Massachusetts; 6, Philadelphia, Pennsylvania; 7, 
Charleston, South Carolina; 8, Bodo, Norway; 9, Bergen, 
Norway; 10, Egersund, Norway; 11, Thurso, Scotland; 12, 
The Hague, The Netherlands; 13, Goteborg, Sweden; 14, 
Karlskrona, Sweden; 15, Visby, Sweden; 16, Helsinki, Fin- 
land. 


Figure 8. As in Figure 4, except that the lines on each 
frame are distinguished by the Earth model used in the cal- 
culations. The solid line was computed using the test Earth 
model (LT = 120 km). The dashed and dotted lines are 
computed by changing LT to 95 and 145 km, respectively. 
The surface load was composed of the ICE-3G adapted ice 
history and a gravitationally self-consistent ocean redistri- 
bution. 


Figure 9. As in Figure 8, except that the dashed and dotted 
lines were computed by changing vlm from the test model 
value (2 x 10 31 Pa s) to 10 31 and 5 x 10 31 Pa s, respectively. 
The arrows labeled A to R indicate the location of 18 specific 
sites on the profiles for which Frechet kernels are presented 
on Figure 12. 


Figure 10. As in Figure 8, except that the dashed and dotted 
lines were computed by changing i/um from the test model 
value (IQ 31 Pa s) to 5 x 10 30 and 3 x 10 30 Pa s, respectively. 
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Figure 11. (top) The present-day Ug component of tangential 
velocities along the great circle profile from Fort Rupert to 
Richmond, relative to the value obtained at Fort Rupert, 
computed using the results in Figure 9. (bottom) As in top, 
except for the European profile and Lulea. 


Figure 12. Frechet kernels for the components of present- 
day velocities computed using the test Earth model. The 
solid, dashed, and dotted lines refer to the kernels for U T , 
Uj, and Ug, respectively. Each frame provides the kernels 
for a different site (see Figure 9) on either the Fort Ruperts 
Richmond profile (A to I), or the Lulea- Wettzell profile (J 
to R). 


Figure 12. (continued) 


Figure 13. The VLBI baselines included in this study (see 
Table 3). The letters denote specific VLBI sites: (A) Algo- 
nquin Park; (B) Beltsville; (E) Effelsberg; (F) Ft. Davis; 
(G) Gilmore Creek; (H) Haystack; (M) Maryland Point; 
(Ma) Madrid; (Me) Medicina; (N) Green Bank; (No) Noto; 
(O) Onsala; (P) Platteville; (R) Richmond; (W) Westford; 
(We) Wettzell; and (Wh) Whitehorse. The shaded region in 
each frame shows the geographic extent of the North Amer- 
ican and Fennoscandian ice cover during the last glacial 
maximum according to the ICE-3G model [Tushingham and 
Peltier, 1991], 


Figure 14. Frechet kernels for a subset of the North Ameri- 
can baselines and all the European baselines generated using 
the test Earth model. The solid, dashed, and dotted lines 
refer to the kernels for the L, V, and T components, respec- 
tively, of the baseline rates of change. 


Figure 14. (continued) 



Figure 15 . The reduced x 3 for the European VLBI baseline 
rate determinations relative to predictions based on three 
families of Earth models. The solid, dashed, and dotted 
lines join results based on varying the value of i>vm , vlm , 
and LT, respectively, in the test model. 
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Figure 1 . Present-day deformation rates computed using the ice disk models with maximum radius (a) 
and (b) 16° and (c) and (d) 8°. The solid line on each frame is generated using the test Barth model. 
The remaining lines are generated using an Earth model constructed by varying a single parameter of 
the test model: LT=96 km (dashed line); vum = 3 x 10 20 Pa s (dashed-dotted line); and vlm = 5 x 10 21 
Pa s (dotted line). As a consequence of symmetry, the tangential velocity is identically zero at the center 
of the disk location and at the antipole. A positive tangential velocity refers to motion directed away 
from the center of the disk load. 


Figure 2 . As in Figures la and lb, except that the dashed line was computed using the test Earth model 
with an elastically incompressible rheology. A positive tangential velocity refers to motion directed away 
from the center of the disk load. 


Figure 3. Plots of the predicted (bottom) present-day radial and (top) tangential velocity over North 
America and Europe computed using the test Earth model and a surface mass load composed of the 
ICE-3G adapted ice history and a gravitationally self-consistent ocean load. The solid dots (top frames) 
refer to the sites: Gilmore Creek, Alaska (G); Fort Rupert, Quebec (F); Algonquin Park, Ontario (A); 
Richmond, Florida (R); Onsala, Sweden (O); Lulea, Sweden (L); and Wettzell, Germany (W). 


Figure 4. (a) The present-day three-dimensional velocities computed along the profile from Fort Rupert 
to Richmond using the test Earth model and a gravitationally self-consistent ocean mass redistribution. 
The solid and dashed lines were computed using ice loads adapted from, respectively, the ICE-3G and 
ICE-1 deglaciation histories, (b) As in Figure 4a, except for the profile extending from Lulea to Wettzell. 


Figure 5. As in Figure 3, except for an ice load adapted from the ICE-1 deglaciation history. 


Figure 6 . As in Figure 4, except all calculations use an ice load adapted from the ICE-3G deglaciation 
history and the test Earth model. The solid, dashed and dotted lines were generated by assuming a 
gravitationally self-consistent ocean mass redistribution, a eustatic ocean mass redistribution, and no 
ocean load component, respectively. 
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Figure 7. Vector plots of present-day tangential velocities at selected coastal sites in (a) eastern North 
America and (b) Europe, computed using the test Earth model and the ICE-3G adapted ice loading 
history. The vectors whose heads are given by solid triangles, open triangles, or two lines, refer to results 
generated using a gravitationally self-consistent, eustatic, or no ocean load. The numbers refer to the 
sites: 1, Fort George, Quebec; 2, Indian Harbour, Newfoundland; 3, Ille d’Anticosti, Quebec; 4, Quebec 
City, Quebec; 5, Boston, Massachusetts; 6, Philadelphia, Pennsylvania; 7, Charleston, South Carolina; 
8, Bodo, Norway; 9, Bergen, Norway; 10, Egersimd, Norway; 11, Thurso, Scotland; 12, The Hague, The 
Netherlands; 13, Goteborg, Sweden; 14, Karlskrona, Sweden; 15, Visby, Sweden; 16, Helsinki, Finland. 


Figure 8. As in Figure 4, except that the lines on each frame are distinguished by the Earth model used 
in the calculations. The solid line was computed using the test Earth model (LT — 120 km). The dashed 
and dotted lines are computed by changing LT to 95 and 145 km, respectively. The surface load was 
composed of the ICE-3G adapted ice history and a gravitationally self-consistent ocean redistribution. 


Figure 9. As in Figure 8, except that the dashed and dotted lines were computed by changing vlm from 
the test model value (2 x 10 21 Pa s) to 10 21 and 5 x 10 21 Pa s, respectively. The arrows labeled A to 
R indicate the location of 18 specific sites on the profiles for which Frechet kernels are presented on 
Figure 12. 


Figure 10. As in Figure 8, except that the dashed and dotted lines were computed by changing vum from 
the test model value (10 21 Pa s) to 5 x IQ 20 and 3 x 10 20 Pa s, respectively. 


Figure 11. (top) The present-day Ug component of tangential velocities along the great circle profile from 
Fort Rupert to Richmond, relative to the value obtained at Fort Rupert, computed using the results in 
Figure 9. (bottom) As in top, except for the European profile and Lulea. 


Figure 12. Frechet kernels for the components of present-day velocities computed using the test Earth 
model. The solid, dashed, and dotted lines refer to the kernels for U r , Uj, and Ug, respectively. Each 
frame provides the kernels for a different site (see Figure 9) on either the Fort Rupert-Richmond profile 
(A to I), or the LuleSrWettzell profile (J to R). 


Figure 12. (continued) 


31 



Figure 13. The VLBI baselines included in this study (see Table 3). The letters denote specific VLBI sites: 
(A) Algonquin Park; (B) Beltsville; (E) Effelsberg; (F) Ft. Davis; (G) Gilmore Creek; (H) Haystack; 
(M) Maryland Point; (Ma) Madrid; (Me) Medicina; (N) Green Bank; (No) Noto; (O) Onsala; (P) Plat- 
teville; (R) Richmond; (W) Westford; (We) Wettzell; and (Wh) Whitehorse. The shaded region in each 
frame shows the geographic extent of the North American and Fennoscandian ice cover during the last 
glacial maximum according to the ICE-3G model [Tushingham and Peltier, 1991]. 


Figure 14. Frechet kernels for a subset of the North American baselines and all the European baselines 
generated using the test Earth model. The solid, dashed, and dotted lines refer to the kernels for the L, 
V, and T components, respectively, of the baseline rates of change. 


Figure 14. (continued) 


Figure 15. The reduced \ 2 for the European VLBI baseline rate determinations relative to predictions 
based on three families of Earth models. The solid, dashed, and dotted lines join results based on varying 
the value of wm, vlm, and LT, respectively, in the test model. 
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Figure 1. Present-day deformation rates computed using the ice disk models with maximum radius (a) and (b) 16° and (c) 
and (d) 8°. The solid line on each frame is generated using the test Earth model. The remaining lines are generated using 
an Earth model constructed by varying a single parameter of the test model: LT=96 km (dashed line); uum = 3 x 10 JO Pa s 
(dashed-dotted line); and ulm = 5 x 10 21 Pa s (dotted line). As a consequence of symmetry, the tangential velocity is 
identically zero at the center of the disk location and at the antipole. A positive tangential velocity refers to motion 
directed away from the center of the disk load. 


Figure 2. As in Figures la and lb, except that the dashed line was computed using the test Earth model with an elastically 
incompressible rheology. A positive tangential velocity refers to motion directed away from the center of the disk load. 


Figure 3. Plots of the predicted (bottom) present-day radial and (top) tangential velocity over North America and Europe 
computed using the test Earth model and a surface mass load composed of the ICE-3G adapted ice history and a gravita- 
tionally self-consistent ocean load. The solid dots (top frames) refer to the sites: Gilmore Creek, Alaska (G); Fort Rupert, 
Quebec (F); Algonquin Park, Ontario (A); Richmond, Florida (R); Onsala, Sweden (O); Lulea, Sweden (L); and Wettzell, 
Germany (W). 


Figure 4. (a) The present-day three-dimensional velocities computed along the profile from Fort Rupert to Richmond 
using the test Earth model and a gravitationally self-consistent ocean mass redistribution. The solid and dashed lines were 
computed using ice loads adapted from, respectively, the ICE-3G and ICE-1 deglaciation histories, (b) As in Figure 4a, 
except for the profile extending from Lule& to Wettzell. 


Figure 5. As in Figure 3, except for an ice load adapted from the ICE-1 deglaciation history. 


Figure 6. As in Figure 4, except all calculations use an ice load adapted from the ICE-3G deglaciation history and the test 
Earth model. The solid, dashed and dotted lines were generated by assuming a gravitationally self-consistent ocean mass 
redistribution, a eustatic ocean mass redistribution, and no ocean load component, respectively. 


Figure 7. Vector plots of present-day tangential velocities at selected coastal sites in (a) eastern North America and (b) 
Europe, computed using the test Earth model and the ICE-3G adapted ice loading history. The vectors whose heads 
are given by solid triangles, open triangles, or two lines, refer to results generated using a gravitationally self-consistent, 
eustatic, or no ocean load. The numbers refer to the sites: 1, Fort George, Quebec; 2, Indian Harbour, Newfoundland; 3, 
Ille d’Anticosti, Quebec; 4, Quebec City, Quebec; 5, Boston, Massachusetts; 6, Philadelphia, Pennsylvania; 7, Charleston, 
South Carolina; 8, Bodo, Norway; 9, Bergen, Norway; 10, Egersund, Norway; 11, Thurso, Scotland; 12, The Hague, The 
Netherlands; 13, Goteborg, Sweden; 14, Karlskrona, Sweden; 15, Visby, Sweden; 16, Helsinki, Finland. 
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Figure 8. As in Figure 4, except that the lines on each frame are distinguished by the Earth model used in the calculations. 
The solid line was computed using the test Earth model (LT = 120 km). The dashed and dotted lines are computed by 
changing LT to 95 and 145 km, respectively. The surface load was composed of the ICB-3G adapted ice history and a 
gravitationally self-consistent ocean redistribution. 


Figure 9. As in Figure 8, except that the dashed and dotted lines were computed by changing vlm from the test model 
value (2 X 10 31 Pa s) to 10 31 and 5 x 10 31 Pa s, respectively. The arrows labeled A to R indicate the location of 18 specific 
sites on the profiles for which Frechet kernels are presented on Figure 12. 


Figure 10. As in Figure 8, except that the dashed and dotted lines were computed by changing uum from the test model 
value (10 31 Pa s) to 5 x 10 30 and 3 x 10 30 Pa s, respectively. 


Figure II. (top) The present-day Ug component of tangential velocities along the great circle profile from Fort Rupert to 
Richmond, relative to the value obtained at Fort Rupert, computed using the results in Figure 9. (bottom) As in top, 
except for the European profile and Lulea. 


Figure 12. Frechet kernels for the components of present-day velocities computed using the test Earth model. The solid, 
dashed, and dotted lines refer to the kernels for U r , U^ and Ug, respectively. Each frame provides the kernels for a different 
site (see Figure 9) on either the Fort Rupert-Richmond profile (A to I), or the Lulea- Wettzell profile (J to R). 


Figure 12. (continued) 


Figure 13. The VLBI baselines included in this study (see Table 3). The letters denote specific VLBI sites: (A) Algonquin 
Park; (B) Beltsville; (E) Effelsberg; (F) Ft. Davis; (G) Gilmore Creek; (H) Haystack; (M) Maryland Point; (Ma) Madrid; 
(Me) Medicina; (N) Green Bank; (No) Noto; (O) Onsala; (P) Platteville; (R) Richmond; (W) Westford; (We) Wettzell; and 
(Wh) Whitehorse. The shaded region in each frame shows the geographic extent of the North American and Fennoscandian 
ice cover during the last glacial maximum according to the ICE-3G model [Tushingham and Peltier, 1991]. 


Figure 14. FYechet kernels for a subset of the North American baselines and all the European baselines generated using the 
test Earth model. The solid, dashed, and dotted lines refer to the kernels for the L, V, and T components, respectively, of 
the baseline rates of change. 
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Figure 14. (continued) 


Figure IS. The reduced x 2 for the European VLBI baseline rate determinations relative to predictions based on three 
families of Earth models. The solid, dashed, and dotted lines join results based on varying the value of vum, »lm, and 
LT, respectively, in the test model. 
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Table 1 . Models Used in Numerical Calculations 


Model 

Lithospheric 

Thickness, 

km 

Upper Mantle 
Viscosity, 
10 21 Pa s 

Lower Mantle 
Viscosity, 
10 21 Pa s 

Test 

120 

1.0 

2.0 

1 

120 

0.3 

2.0 

2 

120 

0.5 

2.0 

3 

120 

1.0 

0.5 

4 

120 

1.0 

1.0 

5 

120 

1.0 

5.0 

6 

71 

1.0 

2.0 

7 

96 

1.0 

2.0 

8 

145 

1.0 

2.0 


All models have an elastic structure given by the seismically 
determined model PREM \Dziewonski and Anderson, 1981). 
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Table 2a. Contributions to Baseline Velocity Predictions of Radial and Tangential 
Site Velocities 




Test Model 


Model 1 



Model 5 


L/T/V 

TOT 

RAD 

TAN 

TOT RAD 

TAN 

TOT 

RAD 

TAN 

L 

r.2 

1.1 

Algonquin Park-Gibnore Creek 
0.2 3.3 0.7 2.7 

0.9 

1.9 

-1.0 

T 

0.7 

0.0 

0.7 

0.5 0.0 

0.5 

1.3 

0.0 

1.3 

V 

4.7 

5.1 

-0.4 

2.0 2.3 

-0.3 

8.3 

8.9 

-0.6 

L 

-1.1 

0.1 

-1.2 

Onsala-WettzeU 
-0.3 -0.1 

-0.4 

-1.3 

0.1 

-1.4 

T 

-0.4 

0.0 

-0.4 

-0.1 0.0 

-0.1 

-0.8 

0.0 

-0.8 

V 

1.3 

1.3 

0.0 

1.0 0.9 

0.1 

2.1 

2.2 

-0.1 


All rates expressed in millimeters per year. L, T, and V indicate the length, transverse, 
and vertical baseline rate vector components, respectively. TOT, RAD, and TAN indicate 
contributions from the total, radial, and tangential site velocities, respectively. 
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Table 2b. Predicted Site Velocities 


Site 

Test Model 

Model 1 

Model 5 


Ur 

Ue 

0 ^ 

U r 

Ue 

Uj, 

Ur 

Ue 

U 4 

Algonquin Park 

4.2 

0.8 

0.1 

2.2 

1.9 

0.4 

7.5 

0.3 

0.1 

Gilmore Creek 

-1.2 

-0.4 

0.4 

-0.2 

0.5 


-2.0 

-1.4 

1.2 

Onsala 

1.2 

0.9 

-0.5 


1.1 

-0.1 

1.7 


-1.1 

Wettzell 

-0.1 

-0.3 

-0.1 


0.7 

0.0 

-0.4 

-1.4 

-0.4 


All rates expressed in millimeters per year. 
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Table 3. Velocities for North American and European Baselines 


Model 


Baseline 


1 

2 

T 

3 

4 

T 

5 

6 

7 

T 

8 

Baseline Rate(s) 

Algonquin Park- 

L 

3.3 

3.0 

1.2 

0.4 

1.0 

1.2 

0.9 

-0.9 

0.4 

1.2 

1.7 

3.3 ±0.8 


Gilmore Creek 

T 

0.5 

0.8 

0.7 

-0.3 

0.0 

0.7 

1.3 

0.8 

0.7 

0.7 

0.7 



— • 

V 

2.0 

2.8 

4.7 

0.7 

1.8 

4.7 

8.3 

4.4 

4.6 

4.7 

4.8 



Algonquin Park- 

L 

0.3 

0.1 

-0.1 

-0.1 

0.0 

-0.1 

-0.1 

0.0 

0.0 

-0.1 

-0.1 

-0.7 ±0.5 

(b) 

Westford 

T 

-0.5 

-0.7 

-0.6 

-0.2 

-0.4 

-0.6 

-0.8 

-0.4 

-0.5 

-0.6 

-0.8 




V 

2.6 

4.0 

5.8 

2.4 

3.7 

5.8 

7.8 

5.2 

5.5 

5.8 

6.0 



Beltsville- 

L 

0.2 

0.1 

-0.1 

0.0 

0.0 

-0.1 

-0.1 

-0.2 

-0.1 

-0.1 

0.0 

3.4 ±1.3 

(a) 

Green Bank 

T 

0.1 

0.1 

0.1 

0.0 

0.1 

0.1 

0.1 

0.1 

0.1 

0.1 

0.1 




V 

-0.1 

-0.2 

-0.2 

-0.3 

-0.3 

-0.2 

0.0 

0.0 

-0.1 

-0.2 

-0.2 



Beltsville- 

L 

-0.7 

-0.8 

-0.5 

0.4 

0.2 

-0.5 

-1.2 

-0.4 

-0.5 

-0.5 

-0.5 

6.9 ± 2.5 


Richmond 

T 

0.4 

0.4 

0.1 

-0.1 

0.0 

0.1 

0.1 

-0.3 

-0.1 

0.1 

0.1 




V 

0.1 

-0.7 

-1.9 

-1.1 

-1.7 

-1.9 

-1.0 

-2.1 

-2.0 

-1.9 

-1.7 



Beltsville- 

L 

0.1 

0.0 

-0.2 

-0.1 

-0.1 

-0.2 

-0.2 

-0.3 

-0.3 

-0.2 

-0.1 

0.0 ± 3.0 

(b) 

Westford 

T 

0.4 

0.4 

0.4 

0.0 

0.1 

0.4 

0.5 

0.3 

0.4 

0.4 

0.4 




V 

-0.3 

-0.2 

-0.5 

0.6 

0.4 

-0.5 

-1.4 

-0.6 

-0.6 

-0.5 

-0.4 



Ft. Davis- 

L 

1.2 

0.6 

-0.8 

0.2 

0.3 

-0.8 

-2.8 

-2.0 

-1.3 

-0.8 

-0.5 

2.7 ± 1.7 


Gilmore Creek 

T 

0.4 

0.3 

-0.2 

0.0 

0.0 

-0.2 

-0.4 

-0.4 

-0.3 

-0.2 

-0.1 




V 

-0.2 

-0.7 

-1.6 

0.3 

-0.6 

-1.6 

-1.8 

-1.6 

-1.6 

-1.6 

-1.5 



Ft. Davis- 

L 

0.4 

0.0 

-0.6 

0.2 

0.2 

-0.6 

-1.7 

-1.3 

-0.9 

-0.6 

-0.4 

1 

H-* 

o 

H- 

o 

(c) 

Green Bank 

T 

0.7 

0.6 

0.0 

-0.4 

-0.3 

0.0 

0.3 

-0.2 

-0.1 

0.0 

0.1 




V 

0.3 

1.0 

2.3 

0.6 

1.6 

2.3 

1.5 

2.5 

2.4 

2.3 

2.1 



Ft. Davis- 

L 

0.8 

0.4 

-0.4 

0.0 

0.1 

-0.4 

-1.4 

-1.4 

-0.8 

-0.4 

-0.2 

-1.1 ±0.3 

(b) 

Westford 

T 

1.1 

1.1 

0.5 

-0.3 

-0.1 

0.5 

0.8 

0.0 

0.3 

0.5 

0.5 

0.9 ±0.9 

(d) 


V 

0.2 

1.2 

2.1 

1.5 

2.3 

2.1 

0.3 

2.0 

2.0 

2.1 

2.0 


- 

Gilmore Creek- 

L 

1.9 

0.8 

-1.8 

-0.3 

-0.5 

-1.8 

-3.6 

-3.6 

-2.5 

-1.8 

-1.4 

-0.5 ± 1.4 

(a) 

Green Bank 

T 

0.1 

0.2 

0.2 

-0.4 

-0.3 

0.2 

0.9 

0.4 

0.3 

0.2 

0.2 




V 

0.3 

0.5 

0.6 

0.7 

0.8 

0.6 

-0.2 

0.8 

0.7 

0.6 

0.5 



Gilmore Creek- 

L 

1.2 

0.5 

-1.2 

0.3 

0.1 

-1.2 

-2.7 

-2.3 

-1.6 

-1.2 

-1.0 

5.6 ± 2.3 


Platteville 

T 

0.0 

0.0 

0.0 

-0.2 

-0.1 

0.0 

0.4 

0.1 

0.1 

0.0 

0.0 




V 

0.5 

0.3 

-0.1 

-■0.4 

-0.5 

-0.1 

-0.1 

-0.2 

-0.3 

-0.1 

0.0 



Gilmore Creek- 

L 

1.2 

0.4 

-1.3 

0.5 

0.4 

-1.3 

-3.9 

-2.7 

-1.9 

-1.3 

-1.0 

1.1 ±1.2 


Richmond 

T 

-0.4 

-0.2 

0.2 

-0.1 

-0.1 

0.2 

0.7 

0.6 

0.3 

0.2 

0.2 




V 

0.1 

-0.3 

-0.9 

0.2 

-0.3 

-0.9 

-1.1 

-0.8 

-0.9 

-0.9 

-0.8 

- 


Gilmore Creek- 

L 

2.4 

1.3 

-1.3 

-0.6 

-0.5 

-1.3 

-2.6 

-3.1 

-2.0 

-1.3 

-0.9 

-0.4 ± 0.5 

(b) 

Westford 

T 

0.1 

0.1 

0.1 

-0.4 

-0.3 

0.1 

0.6 

0.3 

0.2 

0.1 

0.1 

-2.3 ±2.2 

(d) 


V 

0.3 

0.8 

0.6 

1.6 

1.5 

0.6 

-1.1 

0.4 

0.5 

0.6 

0.6 



Gilmore Creek- 

L 

0.1 

-0.1 

-0.5 

-0.1 

-0.1 

-0.5 

-0.9 

-0.7 

-0.6 

-0.5 

-0.4 

-3.8 ± 3.8 


Whitehorse 

T 

0.3 

0.4 

0.2 

-0.2 

-0.1 

0.2 

0.5 

-0.2 

0.1 

0.2 

0.3 




V 

-1.1 

-2.2 

-4.3 

-2.4 

-3.0 

-4.3 

-6.1 

-6.5 

-5.5 

-4.3 

-3.5 
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The model “T” refers to the test Earth model (see text). All rates expressed in millimeters per year. Letters in last column 
are used to indicate which of two antennas at a site (or sites) was involved in obtaining the results shown, according to the 
following scheme: (a) NRAO 85’ antenna, Green Bank; (b) Westford antenna, Westford; (c) NRAO 140’ antenna, Green Bank; 
(d) Haystack antenna, Westford. 
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Table 3. Velocities for North American and European Baselines (continued) 


Model 

Baseline 1 2 T 3 4 T 5 6 7 T 8 Baseline Rate(s) 


o 


41 






&n Jn 

• • 

( jA/uim ) uotjBinjojaQ Xbq juasajj 








eformation Rates ( mm/yr ) 





- 4.5 - 3.0 - 1.5 0.0 1.5 3.0 4.5 6.0 7.5 9.0 10.5 12.0 

mm/yr 


3 















- 4.5 - 3.0 - 1.5 0.0 1.5 3.0 4.5 6.0 7.5 9.0 10.5 12.0 

mm/yr 


S' 















-01 






























h 









120 - 


































3 



l + 20 

10 20 Pas 

Lit! 

1 

1 1 

3xl0 20 

i 

—I 1 1 L 

1 . 

- L - L S 

10 Pas 

t i i 

t- 

t t i 

xio 20 

I0 21 


3xl0 21 


I0 2 

1 

1 

i 




I 


^um(Pos) 


\>lm^ P qs) 

22 Lm 


50 


100 


150 


200 


250 


LT (km) 







